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CONTEUDO 



Prefácio para a IX Escola de 
Computagáo 



Propomos, para a IX Escola de Computagáo, o curso acompanhado de livro texto: Esparsidade, 
Estrutura, Escalamento e Estabilidade em Álgebra Linear Computacional. Embora 
abrangente o suficiente para a compreensáo dos principais problemas da área, a tónica do curso é 
a solugáo de sistemas lineares esparsos e-ou estruturados. Na primeira segáo da introdugáo damos 
uma visáo panorámica sobre o tema e a organizagáo do livro. Neste prefácio ressaltamos algumas 
motivagoes para a inclusáo do curso na IX Escola de Computagáo, e esclarecemos a forma e o 
conteúdo das palestras a serem dadas durante o mesmo. 



0.1 Importáncia da Area 

Grande parte do processamento de dados em ciéncia envolve computagáo numérica, e a maior 
parte desta computagáo numérica sáo rotinas básicas de álgebra linear. Várias aplicagóes em 
engenharia, física, pesquisa operacional, administragáo ou matemática, geram problemas lineares 
cujas matrizes sáo esparsas (entre tantos outros: resolugáo de equagóes diferenciais, análise de 
sistemas, circuitos ou estruturas, programagáo linear e náo linear, otimizagáo de fiuxos em re- 
des, etc...). Ademais, observa-se que a densidade destas matrizes decresce com a dimensáo do 
problema: tipicamente apenas na ou nlog(n) de seus n^ elementos sáo náo nulos. Assim, quanto 
maiores e mais complexas os problemas (e usuários sempre querem resolver modelos maiores), mais 
importante se torna usar eficientemente a esparsidade e a estrutura das matrizes envolvidas. Por 
exemplo, na área de otimizagáo, problemas hoje considerados de grande porte envolvem milhóes 
de equagoes, sendo a maioria destes problemas altamente estruturados, muito esparsos, e usando 
mais de 99% do tempo de resolugáo em rotinas básicas de álgebra linear! 



0.2 Interdisciplinaridade 

Um aspecto que, a nosso ver, torna o tema interessante para um evento como a Escola de 
Computagáo é sua interdisciplinaridade, vejamos: 
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A resolugáo dos problema de álgebra linear envolve os aspectos clássicos de complexidade, 
convergéncia e estabilidade de análise numérica. 

O tratamento de esparsidade e estrutura é um problema essencialmente combinatório, envol- 
vendo teoria de grafos, hipergrafos, e heurísticas para solugáo de problemas de programagáo 
matemática discreta. 

A paralehzagáo destes algoritmos, ou o desenvolvimento de novos métodos, em ambientes 
táo diversos como máquinas de memória compartilhada ou redes de estagóes de trabalho, 
vem hderando as pesquisas na área nos úhimos anos. 



0.3 Serventia do Livro Texto 

O curso que se segue foi recentemente montado como a disciphna MAC-795, Métodos Computa-cio- 
nais da Álgebra Linear, no programa de mestrado do Departamento de Ciéncia da Computagáo 
da Universidade de Sáo Paulo. Cursos de métodos computacionais da álgebra hnear tem se 
popularizado nos úhimos anos em muitos departamentos de computagáo, engenharia e matemática 
aphcada. Era minha intengáo escrever mais um capítulo sobre métodos iterativos mas, ao perceber 
que meu rascunho duphcava o tamanho do presente hvro, resolvi postergar a tarefa para outra 
ocasiáo. 



0.4 Plano de Aulas 

Em cinco palestras de 90 minutos é possível dar um bom panorama da área, sua importáncia, 
problemas, métodos, e perspectivas. O material das palestras foi distribuído da seguinte forma: 

1. Introdugáo: Visáo geral da área, sua importáncia, origem de problemas de grande porte, 
e aspectos essenciais para a sua solugáo. 

2. Esparsidade, EliminaQao Assimétrica: Minimizagáo de preenchimento local, outras 
heurísticas, atuahzagóes de base. 

3. Esparsidade, Eliminagáo Simétrica: Arvores de ehminagáo, heurísticas de ordenagáo, 
heurísticas de dissecgáo. 

4. Estrutura: Esparsidade macroscópica, métodos de blocos, heurísticas de redugáo a formas 
blocadas. 

5. Paralelismo: Uso de esparsidade x uso de estrutura, granularidade, e potenciahdades de 
ambientes com diferentes coeficientes entre velocidade de processamento e velocidade de 
comunicagáo. 
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0.5 Comentário sobre a Bibliografia de Suporte 

Para a parte clássica de análise numérica há uma farta variedade de livros texto, como: [Stewart-73] 
e [Golub-83]. Existem alguns livros ou coletáneas sobre matrizes esparsas em ambiente seqüencial, 
sendo os principais: [Rose-72], [Tewarson-73], [Bunch-76], [Duff-79], [George-81] e também [Pissan- 
84] e [Duff-86]. Com excegáo dos dois últimos, estes hvros já estáo bastante desatuahzados e 
esgotados. 

Álgebra Linear Computacional, sob o enfoque de Computagáo Paralela, é um campo de in- 
tensa pesquisa atuah Exclusivamente para matrizes densas há varias obras pubhcadas, xlcomo: 
[Bertsekas-89] e [Dongarra-91]. Coletáneas de artigos sobre álgebra hnear computacional em am- 
biente paralelo, trazendo alguns artigos sobre matrizes esparsas, sáo muito poucos; como: [Carey- 
89], [Vorst-89] e [Gahivan-90] . Nestas resenhas encontra-se também uma extensiva bibhografia 
comentada da área. 
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Capítulo 1 
INTRODUgÁO 



1.1 Panorama do Livro e seu Tema 

O objetivo deste curso é o estudo de métodos computacionais para a resolugáo de sistemas lineares, 
Ax = b. Enfatizaremos quatro aspectos da construcáo de bons métodos computacionais: 

• Esparsidade: Como resolver eficientemente sistemas cujas matrizes tenham baixa densi- 
dade de elementos náo nulos. 

• Estrutura: Como resolver eficientemente sistemas esparsos cujos elementos náo nulos estáo 
dispostos com uma certa regularidade na matriz de coeficientes. 

• Escalamento: Como minimizar erros de arredondamento gerados ao operarmos com números 
de ordem de grandeza muito diferente, numa representagáo de precisáo limitada. 

• Estabilidade: Como lidar com o efeito cumulativo dos erros de arredondamento na solugáo 
final gerada pelo algoritmo. 

As ferramentas relevantes para a análise e implementagáo de algoritmos eficientes de álgebra 
linear estáo espalhados entre diversas áreas, áreas estas afins-mas-nem-tanto, como Teoria de 
Grafos e Hipergrafos, Álgebra Linear, Anáhse Numérica, e Teoria de Processamento Paralelo. O 
propósito destas notas é apresentar este material de forma razoavelmente coerente e didática. 

O Capítulo 2 exp5e o método de Gauss para solugáo de sistemas hneares, inversáo, ou fatoragáo 
de matrizes. Os métodos para matrizes esparsas hdam com a estrutura da disposigáo dos elementos 
náo nulos dentro da matriz, e esta estrutura é convenientemente descrita e manipulada em termos 
da teoria de grafos. O Capítulo 3 é um resumo de conceitos básicos desta teoria. Conceitos mais 
avangados (ou menos usuais) como grafos cordais, separadores, e hipergrafos, sáo tratados em 
outros capítulos. 

11 
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No Capítulo 4 estudamos técnicas para tratar sistemas esparsos assimétricos, principalmente 
preenchimento local, e a heurística P3. Desenvolvimentos posteriores deste tema encontram-se 
também no capítulo 7. O Capítulo 5 expoem as Fatoragóes QR e de Cholesky, e sua utihzagáo na 
solugáo de problemas de quadrados mínimos, programagáo quadrática, e construgáo de projetores. 
No Capítulo 6 estudamos técnicas para tratar sistemas esparsos simétricos; incluindo toda a 
necessária teoria de grafos cordais. 

O capitulo 7 trata da estrutura de um sistema, que pode ser vista como regularidades no padráo 
de esparsidade, ou como a decomposigáo do sistema em sub-sistemas acoplados. Neste capítulo 
estudamos brevemente a paralehzagáo de alguns dos algoritmos, tema que já aparece imphcita- 
mente em capítulos anteriores. O Capitulo 8 é dedicado á anáhse do efeito dos inevitáveis erros 
de arredondamento nos procedimentos computacionais. No Capitulo 9 anahsamos quanto estes 
erros podem degradar a solugáo ñnal de um problema. Estes dois capítulos tratam dos aspectos 
de Anáhse Numérica que, embora náo sendo a tonica do curso, náo podem ser desprezados. 

O Capítulo 10 trata do problema de mudanga de base, i.é., de atuahzar a inversa de uma 
matriz quando nesta se substitui uma única coluna. Este problema é de fundamental importáncia 
para muitos algoritmos de Otimizagáo. 

Parte da avahagáo numa disciphna como a proposta deve ser feita com exercícios-programa, 
como os dados ao longo das notas. É recomendável o uso de uma hnguagem estruturada, que 
inclua entre seus tipos básicos, números reais de precisáo simples e dupla, inteiros, campos de 
bits e ponteiros, que permita a fácil construgáo e manipulagáo de estruturas, e para a qual haja 
compiladores em computadores dos mais diversos portes. As hnguagens C, C-I--I- e FORTRAN-90 
sáo sugestóes naturais. Ambientes iterativos para cálculo matricial, como Matlab, sáo um grande 
estímulo á experimentagáo e comparagáo de diversos de métodos numéricos, facihtando a rápida 
prototipagem e teste de algoritmos. Apresentamos uma introdugáo a este tipo de ambiente como 
apéndice. 



1.2 Notacoes 

_> 

Utihzaremos letras maiúsculas, A, 5, . . ., para denotar matrizes, e letras minúsculas, a,b, . . ., para 
vetores. Numa matriz ou vetor, os índices de hnha ou coluna seráo, respectivamente, subscritos 
ou superscritos á direita. Assim: Aj é o elemento na i-ésima hnha e na j-ésima coluna da matriz 
A; bi é o elemento na i-ésima hnha do vetor-coluna b; e c^ é o elemento na j-ésima coluna do 
vetor-hnha c. 

índices ou sinais á esquerda, superscritos ou subscritos, identificam vetores ou matrizes distin- 
tas. Assim:A, 'A, ^A, ^A, jA, . . . sáo matrizes distintas. Também a letra S forma, á esquerda de 
uma letra latina o nome de um vetor ou matriz, como 6a ou 6A. As demais letras gregas usaremos 
geralmente para escalares ou fungoes. 

Freqüentemente nos referimos aos blocos componentes de vetores ou matrizes, como por ex- 
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emplo: se 6 e c sao vetores linha a = b c é um vetor linha cujos primeiros elementos sao os 
elementos de b, e os elementos seguintes sáo os elementos de c. 

Analogamente, 



A 



representa uma matriz blocada, desde que as dimensoes dos blocos sejam compatíveis. 

Se A é uma matriz, A^ representa a i-ésima hnha de A, e A^ sua j-ésima coluna, de modo que 
se A é mxn, 



\A 


\A 


\A- 


\A 


\A 


lA 


\A 


lA 


lA 



A 



A^ A^ 



A"-' 



A, 



Ar. 



Uma matriz diagonal será denotada como 

di 
diag{d) = D 



do 



di 
d^ 

(Irt 



O vetor unitário, 1, é o vetor, hnha ou coluna, em que todos os elementos sao iguais a 1, 



isto é, 1 



1 1 



1 



Assim a matriz identidade é I = diag{l). O j-ésimo versor de 
dimensáo n é P , a, j-ésima coluna da matriz identidade de dimensáo n. A transposta, a inversa e 
a transposta da inversa de uma matriz A sáo denotadas, respectivamente, por A' , A~^ e A~*. 

O determinante de uma matriz A, nxn, é 



det{A) = Y.S{p)Al^^Al^,y..A 



p{n) 



onde p 



1 2 

P p 



P" 



é uma permutagáo dos elementos de ^p 



1 2 



n 



O 



sinal de permutagáo, S{p), é +1 ou -1 conforme o número de trocas de pares de elementos que 
é necessário fazer em p para retornar a °p, seja par ou ímpar. O conjunto dos elementos em 
cada um dos termos na somatória da definigáo do determinante é denominado uma diagonal da 
matriz. A diagonal correspondente a pemutagáo ^p é denominada diagonal principal. Uma 
matriz quadrada é singular se tiver determinante nulo. O posto de uma matriz A, rank{A), é a 
dimensáo de sua maior sub-matriz quadrada náo singular. 

Dado um sistema Ax = b, com matriz de coeficientes A nxn náo singular, a Regra Cramer 
nos diz que o sistema tem por única solugáo 



X 



Xi 



det{ 



A' 



A 



i-l 



b A 



i+l 



^n 



)/det{A). 
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O trago de uma matriz A, nxn, é a soma de seus elementos na diagonal principal: 

n 



i=l 



A matriz booleana associada á matriz A, B{A), é a matriz, da mesma dimensáo de A, em 
que B{A)l = 1 se A^ 7^ 0, e B{A)l = se A^ = 0. O complemento de uma matriz booleana B, é 
a matriz booleana 5, da mesma dimensáo de -B, tal que B¡ = 1 -^ B¡ = 0. 

Os conjuntos A^ = {1, 2, ..., n} e M = {1, 2, ..., m} seráo freqüentemente usados como domínios 
de índices. Assim, se A é uma matriz mxn, faz sentido falar dos elementos Aj, i & M, j G A^, 

O número de elementos náo nulos, ENNs, numa matriz A, mxn, é 

m,n 

enn{A) = ^ B{A)¡ = l'B{A)l 

de modo que enn{Ai) e enn{A^) sáo, respectivamente, o número de elementos náo nulos na i-ésima 
linha e na j-ésima coluna de A. 

Dada uma fungáo (/?(), definida num domínio D, e X G D, definimos seu argumento mínimo 
V = arg minxex <^{x) e argumento máximo U = argmaxxex y^{x) como, respectivamente, os 
conjuntos V, U G D que minimizam ou maximizam a fungáo ^p em X. Assim, se por exemplo, 

(/9 = x^, X G 3?, X = [—5,5] e Y =] — 5,5[, entáo argmmxVÍ.^) = argmm.Y f^^x) = {0}, 
argm?íXx^{x) = {—5,5}, argm?íXY^{x) = 0. 

Para realizar experiéncias computacionais utilizaremos por vezes os sistemas lineares de di- 
mensáo n e solugáo x = 1, com as seguintes matrizes de coeficientes e vetores independentes: 



Binomial: 



Hilbert: 



Tridiagonal: 



nrpj 



""B^ 



l 

J-1 





^ j = l,2,...i 
^ 3>i 



21 
22 

2" 



^HÍ = l/{i+j-r) , -h, = Y,l/{r + j-l] 



T 



-2 ^ 


i = j 


-1 ^ 


\i-j =1 





caso contrario 



1 <í^ i e {l,n} 
caso contrario 
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1.3 Representacáo de Matrizes Esparsas 

Ao representarmos uma matriz esparsa gostaríamos de utilizar o mínimo de memória, idealmente 
apenas a posigáo e o valor dos ENNs, mas ao mesmo tempo ter a maior flexibilidade possível para 
acessar, modiflcar, inserir ou remover um elemento qualquer. Estes objetivos sáo algo conflitantes, 
o que motiva o uso de diversas representagóes: 

Representagáo estática por linhas: Para uma matriz A, mxn com enn{A) = l, sáo utilizados 
um vetor real, aias{k) k E L, um vetor inteiro aijs{k), k E L = {1, . . .1}, e um segundo vetor de 
inteiros, aif{i), i G M. Os vetores aias e aijs listam os valores e os índices de coluna dos ENN's, 
linha por linha. aif{i) nos dá o flm, ou a posigáo do úhimo elemento da hnha i em aias e aijs. 

Representagáo estática por colunas: Para uma matriz A, m x n com enn{A) = l, sáo 
utihzados um vetor real, ajas{k) k E L, um vetor inteiro ajis{k), k E L, e um segundo vetor de 
inteiros, ajf{j), j G A^. Os vetores ajas e ajis listam os valores e os índices de linha dos ENNs, 
coluna por coluna. ajf{i) nos dá a posigáo do úhimo elemento da coluna j em ajas e ajis. 

Representagáo de hsta ligada por hnhas: Cada ENN corresponde a uma célula contendo: O 
valor do ENN, o índice de coluna e um ponteiro para a célula do próximo ENN na hnha. As células 
dos úhimos ENNs em cada linha contém o ponteiro nulo, e um vetor de m ponteiros, ancorai{i) , 
nos dá a primeira célula de cada linha. 

Representagáo de lista hgada por colunas: Cada ENN corresponde a uma célula contendo: 
O valor do ENN, o índice de linha e um ponteiro para a célula do próximo ENN na coluna. As 
células dos úhimos ENNs em cada coluna contém o ponteiro nulo, e um vetor de n ponteiros, 
ancoraj{j), nos dá a primeira célula de cada coluna. 

Representagáo de rede: Cada ENN corresponde a uma célula contendo o valor do ENN Aj, 
os índices i e j, e ponteiros para a célula do próximo ENN na linha e na coluna. Dois vetores 
ancorai{i) e ancoraj{j) contém ponteiros para as primeiras células em cada hnha e coluna. 

Representagáo de rede dupla: Análoga a representagáo de rede, mas cada célula contém além 
de ponteiros para as próximas células ao sul (próximo ENN na coluna) e ao leste (próximo ENN 
na hnha), também ponteiros para as células ao oeste e ao norte, i.e. para os ENN anteriores na 
hnha e na coluna. 

No exemplo 1 temos as diversas representagóes da matriz: 



A 





102 




104 








201 














301 


304 




405 


, / = 


= enn{A) = 


= 9 


501 


502 


503 
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Representagao estática por linhas: 



aias 



102 104 201 301 304 405 501 502 503 



aijs 



2 4 114 5 12 3 



aif 



2 3 5 6 9 



Representagao estática por colunas: 



ajas 



201 301 501 102 502 503 104 304 405 



ajis 



235155134 



ajf 



3 5 6 8 9 



Representagao de lista ligada por linhas: 

ancorai(l) -^ (2, 102, ->) (4, 104, H) 

ancorai(2) -^ (1, 201, H) 

ancorai(3) -^ (1, 301, -^) (4, 304, H) 

ancorai(4) -> (5, 405, H) 

ancorai(5) -^ (1, 501, ^) (2, 502, ^) (3, 503, 

Representagáo de hsta hgada por colunas: 



ancoraj(l) ancoraj(2) ancoraj(3) ancoraj(4) ancoraj(5) 



4- Y 4- 

(2, 201, i) (1, 102, i) (5, 503, 

(3, 301, ;) (5, 502, ±) 

(5, 501, ±) 



Representagao em rede: 
ancoraj(l) 



ancorai(l) — )> 
ancorai(2) — )> 
ancorai(3) — )■ 
ancorai(4) — )■ 
ancorai(5) — )■ 



i 



2 1 201 

i H 

3 1 301 



5 1 501 



ancoraj (2) 

1 2 102 

i ^ 



5 2 502 



i i 

(1, 104, i) (4, 405, 
(3, 304, ±) 



ancoraj (3) 



5 3 503 



ancoraj (4) 

i 
1 4 104 

i ^ 



3 4 304 



ancoraj(5) 
4 



4 5 405 
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Exercícios 

1. Prove que 

(a) {ABy = B'A'. 

(b) (AB)-^ = A-^B-\ 

(c) (A')-^ = {A-^y. 

2. Na representagáo estática por linhas, é realmente necessário termos aif além de aias e aijs7 

3. Considere uma matriz esparsa de estrutura regular e conhecida, por exemplo tri-diagonal, 
para a qual bastaria conhecermos o valor dos ENNs numa dada seqüéncia, por exemplo hnha 
por hnha. Suponha que um inteiro ou ponteiro ocupa 2 bytes e que um real ocupa 4 bytes. 
Dé o coeficiente de uso de memória de cada uma das outras representagóes em relagáo a está 
representagáo minimaL 

4. (a) Considere uma matriz esparsa A de densidade, i.e. fragáo de elementos náo nulos, a"^. 

Suponha que A está representada em rede. Queremos adicionar um novo ENN Aj á 
matriz. Supondo que os ENN estáo aleatoriamente distribuídos, e tomando acesso a 
uma célula como operagáo elementar, qual a complexidade esperada desta operagáo? 
Exphque sucintamente como reahzar a operagáo. 

(b) Nas mesmas condigóes, queremos substituir duas hnhas A^ e Aj+i, respectivamente, 
pelas combinagóes hneares Ai = c* Ai + s * Aj+i e Aj+i = c* Ai^i ~ s* A^. Novamente 
queremos um algoritmo, descrito sumária mas claramente, e a complexidade esperada. 
Dicas e curiosidades: 

i. Assumindo que « << 1, o caso em que estamos interessados, a densidade esperada 

das novas hnhas é 2 * a. Por qué? 
ii. O algoritmo para a segunda questáo deve ser algo melhor que a mera repetigáo do 

algoritmo da primeira questáo. 
iii. Esta transformagáo hnear, tomando as constantes c e s como o coseno e o seno de 

um ángulo 9, é uma "rotagáo de Givens", a ser estudada no capítulo 3. 

5. Dados u e w n X 1, AeBnxn, ek^ N*, prove que 

(a) tr{A + B) = tr{A) + tr{B). 

(b) tr{AB) = tr{BA). 

(c) tr(-uw') = w'u. 

(d) tr{Auw') = tT{uw'A) = w'Au. 

(e) {uw'A)'' = {tT{uw'A))''-^ {uw'A). 

(f) {Auw')'' = {tT{Auw'))''-^ {Auw'). 



18 CAPITULO 1. INTRODUQAO 



Capítulo 2 
FATORAgÁO LU 



Solugáo de Sistemas Lineares 



2.1 Permutacoes e Operagoes Elementares 

Uma matriz de permutaQao é uma matriz obtida pela permutacáo de linlias ou colunas na 
matriz identidade. Realizar, na matriz identidade, uma dada permutagáo de linhas, nos fornece 
a correspondente matriz de permutacáo de linhas; Analogamente, uma permutagáo de colunas da 
identidade fornece a correspondente matriz de permutagáo de colunas. 

Dada uma (matriz de) permutagáo de hnhas, P e uma (matriz de) permutagáo de colunas, Q, 
o correspondente vetor de índices de hnha (coluna) permutados sáo 



p={P 



1 
2 

m 



1 2 



n 



Q 



Lema 2.1 Realizar uma permutagáo de linhas (de colunas) numa matriz qualquer A, de modo a 
obter a matriz permutada A, equivale a multiplicá-la, d esquerda (d direita), pela correspondente 
matriz de permutagdo de linhas (de colunas). Ademais, se p (q) é o correspondente vetor de 
índices de linha (de coluna) permutados. 



M = iPA)¡ 



Kii) 



ÁÍ = {AQ)¡ = A 



lij) 
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Exemplo 1: Dadas as matrizes 



A 



11 


12 


13 ' 


21 


22 


23 


31 


32 


33 



p 



p = q 



3 12 



PA 





1 


= 


1 




1 


31 


32 33 


11 


12 13 


21 


22 23 



Q 



AQ 






1 


' 








1 


1 









13 11 12 
23 21 22 
33 31 32 



Uma matriz quadrada, A, é simétrica sse for igual a transposta, isto é, sse A = A'. Uma 
permutaQao simétrica de uma matriz quadrada A é uma permutagáo da forma A = PAP', 
onde P é uma matriz de permutagáo. Uma matriz quadrada, A, é ortogonal sse sua inversa for 
igual a sua transposta, isto é, sse A~^ = A'. 

Lema 2.2 (a) Matrizes de permutagáo sáo ortogonais. (h) Uma permutagáo simétrica de uma 
matriz simétrica é ainda uma matriz simétrica. 



Estudaremos a seguir alguns algoritmos para solugáo do sistema Ax = b. Tais algoritmos sáo 
denominados diretos pois, a menos do erro de arredondamento, fornecem diretamente a solugáo do 
sistema. A esséncia destes algoritmos sáo transformagóes sobre o sistema que deixam inalteradas 
sua solugáo. 

Sáo operagoes elementares sobre uma matriz, n x m, qualquer: 

1. Multiplicar uma linha por um escalar náo nulo. 

2. Adicionar, a uma linha, uma outra hnha da matriz. 

3. Subtrair de uma hnha, uma outra hnha da matriz muhiphcada por um escalar náo nulo. 

Uma operagáo elementar aphcada á matriz identidade I, n x n, produz a matriz elementar 
correspondente, E. 

Lema 2.3 Realizar uma operagáo elementar sobre uma matriz qualquer, A, equivale a multiplicá- 
la, á esquerda, pela matriz elementar correspondente. 

Exemplo 3. 



A 



EA 



11 12 


13 






1 





21 22 


23 




E = 


1 





31 32 


33 






-31/11 


1 


" 11 


12 




13 




21 




22 




23 





32-12*31/11 33-13*31/11 
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Lema 2.4 Toda matriz elementar é inverswel. Ademais, Ax = b -^ EAx = Eb. 

A matriz aumentada correspondente ao sistema Ax = b, Anxnéa matriz [A b], n x n + 1. 
Obviamente a matriz aumentada do sistema EAx = Eb é a matriz E[A b]. 

Uma matriz quadrada A é dita triangular superior se todos os elementos abaixo da diagonal 
principal sáo nulos, e é dita triangular estritamente superior se também os elementos da diagonal 
sáo nulos. Analogamente, definimos matriz triangular inferior e estritamente inferior. Assim, A é 

• Triangular Superior -^ {i > j ^ Aj = 0). 

• Triangular Estritamente Superior -^ {i > j ^ Aj = 0). 

• Triangular Inferior -^ {i < j ^ Aj = 0). 

• Triangular Estritamente Inferior -^ {i < j ^ Aj = 0). 

Os métodos diretos que estudaremos funcionam por ser fácil resolver um sistema Ux = b se U 
é triangular superior, isto é Em um tal sistema podemos calcular por substituigáo, nesta ordem, 
as componentes x„, x„_i, . . . ,xi, pois 



Xn—k 



bn/K 



{bn~k — 2^ 

j=n—k+l 



^n-k^j)l^n-k 



Tal método funciona se f/j ^ 0, Vj' G A^. Como det{U) = Y\!j=iUj, esta condigáo equivale a 
termos um sistema bem determinado. Estudaremos a seguir como levar, através de operagóes 
elementares, um sistema qualquer á forma triangular. 



2.2 Método de Gauss 

O método de Gauss leva o sistema á forma triangular através de operagóes elementares do tipo 
3. Tentaremos fazé-lo usando, nesta ordem, a primeira linha para anular os elementos abaixo 
da diagonal na primeira coluna, a segunda linha para anular os elementos abaixo da diagonal na 
segunda coluna, etc... O exemplo 2 ilustra o processo. Indicamos também os fatores pelos quais 
multiplicamos cada linha antes de subtraí-la de outra. 

Exemplo 4. 



°A°6 



2 


1 


3 


1 




2 


3 


6 


2 


1 


4 


4 


6 


6 


2 



-)■ 



'Ah 



2 


1 3 


1 " 





2 3 


1 





2 


4 



-)■ 
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'A% 



' 2 1 


3 1 " 




1 


2 


3 1 


-í- X = 


2 





-3 3 




-1 



Em cada etapa denominaremos a linha que está sendo multiplicada e subtraída ás demais de 
linha pivó, seu elemento diagonal de elemento pivó e os fatores de multiplicagáo de multi- 
plicadores. Posteriormente faremos uso dos multiplicadores utilizados no processo e, portanto, 
devemos armazená-los. Com este fim, definimos as matrizes, nxn, ^M, ^M, . . . "^-"^M = M, onde: 
Sej<kei>j, entáo '^M¡ é o multiplicador utilizado na /c-ésima etapa para anular o elemento 
na i-ésima linha e j-ésima coluna; Caso contrario, ^M/ = 0. 

Observe que os elementos náo nulos de ''M correspondem a elementos nulos em '^A, e vice-versa. 
Podemos pois poupar memória, guardando uma única matriz, '^A + ^M. Assim, no Exemplo 4, 
'^A+'^ M I ^6 , k = 0,1,2, pode ser representado como: 



2 1 


3 1 " 




2 3 


6 2 


-)> 


4 4 


6 6 





2 13 1 
i 2 3 1 
^204 



—?■ 



2 1 3 1 
12 3 1 
^ i -3 3 



2.3 Pivoteamento 



Consideremos a hipótese do surgimento de um zero na posigáo do pivo da etapa seguinte do 
processo de triangularizagáo, isto é, em '^^^A anula-se o elemento ''"^Al. Se na coluna '^-^A'^ 
houver algum elemento náo nulo na hnha l, l > k, podemos permutar as hnhas k e l e continuar 
o processo de triangularizagáo. Uma permutagáo de hnhas, para trocar o elemento pivó a ser 
utihzado, denomina-se um pivoteamento. A cada etapa queremos lembrar-nos de quais os 
pivoteamentos realizados durante o processo. Para tanto, guardamos os vetores de índices de 
hnha permutados da permutagáo corrente em relagáo ao sistema originaL Estes sáo os vetores de 
permutagáo, ^p, "^p, ... "~^p = p. No que tange ao armazenamento dos muhiphcadores, devemos 
convencionar se estes seráo ou náo permutados, junto com as respectivas hnhas de A + M, nas 
operagoes de pivoteamento. Adotaremos, por hora, a convengáo de SIM, permutá-los, junto com 
os pivoteamentos. 

O exemplo 5 ilustra o processo para uma matriz de dimensáo k = 4, apresentando a matriz 



'A+'' M 



juntamente com o vetor de permutagao, '^p para cada etapa da triangularizagao: 



2 1 9 


-1 


1 


3 


9 6 


6 


4 


1 3 7 

2 8 4 


7 
2 


2 

3 ^ 


1/3 
2/3 


5 
2 


5 
-2 


2 
3 


3 9 6 


6 . 


4 


. 2/3 


-5 5 


-5 _ 


1 



-)■ 
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3 


9 6 


6 


4 


3 


9 6 


6 


4 


2/3 
2/3 


-5 5 
-2/5 2 


-5 
-4 


1 

3 ^ 


2/3 
1/3 


-5 5 
5 


-5 
5 


1 
2 


1/3 


5 


5 _ 


2 


. 2/3 


-2/5 2/5 


-6 . 


3 



Observagáo 2.1 Note que se o sistema é bem determinado, um pivoteamento que permite o 
prosseguimento do processo de triangularizagáo é sempre possível: Se na /c-ésima etapa, para 
/ = k...n, ^~^Af = 0, isto é o elemento na posicáo pivo e todos os elementos abaixo dele 
na coluna ^~i^^ se anulam, entáo as linhas ^~^Ai, l = k...n, sáo linearmente dependentes e 
det{ ''~^A) = ^ det{ ^A) = 0, o que contraria a hipótese do sistema ser bem determinado. 



Observagáo 2.2 Assuma sabermos de antemáo um vetor de permutagáo, ^~^p = p, viável no 
processo de triangularizagáo de um dado sistema [^A °6]. Neste caso, poderíamos permutar as 
hnhas do sistema original como indicado no vetor de permutagáo, obtendo um novo sistema 
P [^A ^b], cujas equagóes sáo as mesmas que as do sistema original, a menos da ordem em que estáo 
escritas. Poderíamos entáo triangularizar este sistema sem necessidade de nenhum pivoteamento, 
obtendo ao final o mesmo sistema equivalente, ["~^A "^^6]. 

Note que a triangularizagáo da matriz dos coeficientes de um sistema é completamente inde- 
pendente do vetor dos termos independentes. Suponha termos triangularizado o sistema ['^A ^b], 
isto é, que temos um sistema equivalente e triangular ["~^A = U "'~^b], tendo sido preservados os 
multiphcadores e as permutagoes utihzadas, M e p. Se for necessário resolver um segundo sistema 
[^A '^b], que difere do primeiro apenas pelo vetor dos termos independentes, náo será necessário re- 
triangularizar a matriz A; Bastará efetuar no novo vetor de termos independentes, °c, as operagoes 
indicadas pelo vetor de permutagáo p e pela matriz de muhiphcadores M. 

Uma poh'tica, que mais tarde demonstraremos útil, é reahzarmos pivoteamentos para colocar 
como elemento pivo, em cada etapa, o elemento de máximo módulo. Esta estratégia, o pivotea- 
mento parcial, garante termos todos os muhiphcadores | M/ | < 1. O pivoteamento parcial é 
de fundamental importáncia para a estabihdade numérica do processo, controlando a propagagáo 
de erros de arrendondamento. 

Observagáo 2.3 Para reahzar uma operagáo de pivoteamento náo é necessário efetivamente per- 
mutar hnhas da matriz A: Basta, na /c-ésima etapa da triangularizagáo utihzar o índice ^~^p{i) 
ao invés do índice de hnha i. 



Observagáo 2.4 Uma maneira ahernativa de guardar pivoteamentos reahzados ao longo do pro- 
cesso de triangularizagáo é o vetor de pivoteamentos, t, onde t{j) = i significa que ao ehmin- 
armos a coluna j permutamos para a posigáo de pivo (hnha j) o elemento na hnha i. O vetor de 



pivoteamentos do exemplo anterior é 



4 12 3 



. Com o vetor de pivoteamentos adotaremos 



a convengao de NÁO permutar os muhiphcadores na coluna ^ M, nos pivoteamentos k > j. 
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2.4 Lema da Fatoracáo 

_> 

Lema 2.5 (da Fatoragáo) Seja ^A uma matriz inverswel triangularizável pelo método de Gauss 
sem nenhum pivoteamento e sejam, conforme a nomenclatura adotada, A = °A, ^A, . . . ""^A, e 
^M, ^M, . . . "-"^M = M, respectivamente, a k-transformada da matriz A e a k-ésima matriz dos 
multiplicadores. Fazendo U = "'~^A e L = M + I , temos que L e U sáo matrizes triangulares, 
respectivamente inferior e superior, e que A = LU. 

Demonstracáo. 

Verifiquemos inicialmente que a transformacáo linear que leva ^^^A em ^A é dada por ^T, i.e. 
fc^ =fc T ^-M, onde 

" 1 



^rji 



• 



-M^i 



••• 
1 



-M„^ 

de modo que U = ""^/1 = "--^T ^-'^T ... '^T ^T A = TA. Observando ainda que ^T é inversível 
e que 



krp — l k T 



1 
••. 







+M^ 



fc+1 







+M^ 



••• 

1 



Ademais, temos que T ^ = ^T ^ ^T ^ . . . " ^T ^ = ^T . . . " ^L e é fácil verificar que ^L ^L . . . " ^L 
L, donde A = LU. Q.E.D. 

Teorema 2.1 (LU) Seja A uma matriz inversível n x n. Entáo existe uma (matriz de) per- 
mutagáo de linhas P de modo que A = PA é triangularizável pelo método de Gauss e A = LU. 



Demonstracáo. 

O teorema segue trivialmente do lema da fatoracáo e das observagoes 2.1 e 2.2. 

A solucáo de um sistema linear Ax = b pode ser expressa diretamente em termos da fatoracáo 
LU da matriz dos coeficientes, i.e., se A = PA = LU, entáo P'LUx = b, e x = U^^L^^Pb. Assim, 
em muitas aplicagóes, o conhecimento explícito de A~^ pode ser substituído com vantagem pelo 
conhecimento da fatoracáo da matriz. 
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Observagáo 2.5 Considere os dois algoritmos seguintes, o de substituigáo e o de multilpicagáo 
pela inversa, para solucáo do sistema Ux = b: O primeiro algoritmo acessa a matriz U por linha, 
enquanto o segundo acessa a matriz U por coluna. 

X = b ; 

x(n) = x(n) / U(n,n) ; 

for i = n-l:-l: 1 

x(i) = ( x(i) - U(i,i+l:n)*x(i+l:n) ) / U(i,i) ; 
end 

X = b ; 

for j = n:-l:2 

x(j) = x(j) / U(j,j) ; 

x(l:j-l) = x(l:j-l) - x(j)*U(l:j-l,j) ; 
end 
x(l) = x(l) / U(l,l) ; 



2.5 O Método de Doolittle 

Usando o Teorema LU podemos determinar L e U diretamente da equacáo de decomposigáo, 
LU = A, tomando, nesta ordem, para i = 1 . . . n, as equagóes das componentes de A na i-ésima 
linha e na i-ésima coluna, isto é, 

. , { L,U = A, 

para i = l . . .n { 
^ \ LU' = A' 

temos 

. _ í para ]=i...n UÍ = A^ - EÍ=\ M^UÍ 

P''''''" ■■■''l paraj=^ + l...n M^ = {A^ -JX^M^U^/UÍ 

Note que só escrevemos as equagóes para os termos incógnitos, isto é, ou acima da diagonal 
em U, e abaixo da diagonal em L. Também as somatórias foram interrompidas quando os termos 
remanescentes sáo todos nulos. 

O cálculo do elemento (M + f/)¿ , envolve, na ordem prescrita, apenas elementos já calculados 
de (M+f/). Ademais, o cálculo de (M+f/)¿ é a úhima ocasiao em que se necessita os elementos Al; 
portanto podemos armazenar {M + U)l no lugar de Aj. Estas sáo as matrizes A = °D, . . . , "'D = 

M + U. 

Exemplo 6 - Usando o método de Doohttle para triangularizar a matriz do Exemplo 4, onde 
náo há necessidade de pivoteamento, temos 
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''D = A 



2 1 3 
2 3 6 

4 4 6 



-)> 
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6 
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-)> 
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1 3 ' 


i 


3 6 


^ 


4 6 



1/} 



-)■ 



^ 


1 


3 ' 


1 


2 


3 


2 


1 


6 



''D ^ 



2 1 3 

1 2 3 

2 1 -3 



W = M + U 



Para realizar os pivoteamentos necessários á passagem de ''D a '^+^D é necessário examinar os 
possíveis elementos pivos nas linhas i = k . . .n, isto é, os elementos em '^A''. Para tanto, basta 
calcular, para i = k . . .n, 

fc-i 



A'l - Y: M¡Ut . 



1=1 



Exemplo 7 - Triangularizando pelo método de Doolittle a matriz, A, com pivoteamento parcial, 
temos os '^p, ''D, '^v, para k = . . .n, como segue: 
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2.6 Complexidade 

Contemos agora o número de operagóes aritméticas necessárias á triangularizagáo de uma matriz 
de coeficientes e á solugáo de um sistema linear. Na fase de triangularizagáo da matriz dos 
coeficientes, vemos que o cálculo de ^A a partir de '^'^A requer (n — k) divisóes e {n — k)^ somas 
e produtos. No total sáo portanto requeridos 

ji-i 



^(n — /c) = n{n — l)/2 divisoes, e 
fc=i 

^(n — A;)^ = n^n'^ — l)/3 + n{n — l)/2 produtos e subtragóes. 



fc=i 

n-l 



fc=i 
isto é, sáo necessários da ordem de n^/3 + 0{n'^) produtos e subtragóes, e n^/2 + 0{n) divisoes. 

Analogamente, dada a matriz dos multiplicadores e o vetor das permutagóes, M e p, o trata- 
mento de um vetor de termos independentes requer n(n — 1)/2 produtos e subtragóes. Finalmente, 
a solugáo do sistema linear triangularizado requer n divisoes e n{n — l)/2 produtos e subtragóes. 

Exercícios 

1. Comente como cada forma de representagáo da matriz A, por linhas ou por colunas, favorece 
o uso de um dos algoritmos apresentados na observagáo 2.5. Escreva algoritmos similares 
para a solugáo do sistema Lx = b. 

2. Programe e implemente, em linguagem C, C++, ou FORTRAN-90, fungóes para: 

(a) Fatoragáo LU com pivoteamento parcial. 

(b) Solugáo dos sistemas Lx = b e Ux = b, 

(c) Um argumento adicional deve indicar a forma de representagáo das matrizes, densa, 
estática por linha, ou estática por coluna. 

(d) No caso da representagáo densa, um segundo argumento deve indicar o uso do vetor 
de permutagoes ou pivoteamentos. No caso das representagóes estáticas, escreva o 
programa como Ihe parecer mais conveniente. 

3. Mostre que calcular exphcitamente a inversa de A imphca saber resolver n sistemas hneares, 
Ax = P. Dada a fatoragáo A = LU , qual o custo de computar exphcitamente a inversa 
A~^l 
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Capítulo 3 



RESUMO DE TEORIA DOS GRAFOS 



3.1 Conceitos Básicos 



Um grafo é um par ordenado G = [Vg, Tg) onde o primeiro elemento é um conjunto finito, o 
conjunto de vértices, e o segundo elemento é uma fungáo Tg '■ Vg ^-> P(Vg'), no conjunto das partes 
de Vg, a fungáo de filiaQao. Tomaremos Vg um conjunto indexado, Vg = {f i, f^, • . . fn} ou entáo 
tomaremos Vg como sendo o próprio conjunto dos n primeiros inteiros positivos, N = {1,2, . . .n}. 
Para náo sobrecarregar a notagáo escreveremos, quando náo houver ambigüidade, {Vg, Tg) como 

É comum representarmos um grafo por conjunto de pontos (os vértices) e um conjunto de 
arestas (setas) que váo de cada vértice para os seus filhos, isto é, de cada v & V para os elementos 
em r(f ). Podemos definir o grafo através dos seus vértices e de suas arestas, G = {Vg, Ag), onde 
cada aresta é um par ordenado de vértices e {i,j) G Ag ^ j G Pg^í). Uma terceira maneira de 
definir um grafo é pelo par G = {Vg, Bg) onde Bg é a matriz de adjacéncia, a matriz Booleana 
tal que B¡ = l ^ j e V{i). 

Exemplo 1: 

Considere o grafo de vértices A^ = {1,2,3,4,5,6} e fungáo de fihagáo r(l) = 0, r(2) = {2}, 
r(3) = 0, r(4) = {3,5,6}, r(5) = {3,4,5} e r(6) = {5}. Suas arestas e matriz de adjacéncia sáo, 
Ag = {(2, 2), (4, 3), (4, 5), (4, 6), (5, 3), (5, 4), (5, 5), (6, 5)}, e 



5, 



G 





10 



10 11 

1110 

10 



1 3 ^ 5^ 

t /v/ t 



^2 



4 



-)■ 



6 



Dado um subconjunto de vértices W C V, definimos r(W^) = Uw(z]vT{w). Definimos também 
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W, rM 



r(i) e, para A; > 1, F* 



r(r ^(i)). Estes sao, para k = 1,2,3,... 



os filhos, netos, bisnetos, etc. do vértice i. Finalmente definimos os descendentes de i por 
f (i) = U^^qT''{í), e o fecho transitivo G = {Vg,Tg). 

Exemplo 2: 

Damos a matriz de adjacéncia de um grafo, e do respectivo fecho transitivo: 
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Dado o grafo G = {V,r) definimos a fungáo de paternidade i E T ^{j) -^ j E T{i). 
Definimos também seu grafo inverso, G^^ = {V,T^^). 

Lema 3.1 Dado um grafo definido por sua matriz de adjacéncia, G = {N,B), seu grafo inverso 
éG-^ = {N,B'). 



Um sub-grafo de G = {V, A) é um grafo G' = {V', A'), onde V' G V e A' é um subconjunto 
de arestas de A que tem ambos os vértices em V'. O sub-grafo induzido por um subconjunto 
de vértices V' é G = {V', A') onde A' é máximo, e o sub-grafo induzido por um subconjunto de 
arestas A' é G = {V', A') onde V' é mínimo. 

Uma aresta que parte e chega no mesmo vértice é dita um loop. Duas arestas, {i,j) e {k, h), 
sáo ditas contíguas se a primeira chega no vértice de que parte a segunda, isto é se j = k. Um 
passeio é uma seqüéncia náo vazia e finita de arestas contíguas. Um circuito é um passeio que 
parte e chega no mesmo vértice, isto é, o primeiro vértice da primeira aresta é o segundo vértice 
da úhima aresta. Uma trilha é um passeio onde náo se repete nenhuma aresta e um caminho 
é uma trilha na qual náo há (subseqüéncias que sejam) circuitos. Uma trilha na qual o único 
circuito é toda a trilha é uma ciclo. 

Notemos que um passeio C = {wq, Wi), {wi, w^), . . . {w^-i, w^)), é igualmente bem determinado 
pela seqüéncia de seus vértices C = {wq, Wi, w^, . . . w^). 

Teríamos assim, no Exemplo 1, exemplos de: 

. loop: ((2,2)). 

. passeio: ((4, 5), (5, 5), (5, 5), (5,4), (4, 5). 

. circuito: (5,5,4,6,5,5). 

. trilha: (5,4,6,5,5,3). 

. caminho: (6,5,4,3). 
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• ciclo: (5,4, 6, 5)o-u(5, 5). 

Um grafo é acíclico se náo contém ciclos. Uma árvore de raiz v é um grafo acíclico, H = 
{Vh, ^h) , V G Vh, onde todos os vértices tém no máximo um pai e apenas a raiz náo tem pai. Os 
vértices sem filhos de uma árvore dizem-se folhas da árvore. Uma colegáo de árvores é denominada 
floresta. Uma fioresta cobre um grafo G sse é um subgrafo de G que contém todos os seus vértices. 
Seguem exemplos de fiorestas que cobrem o grafo do Exemplo 1, estando assinaladas as raízes. 

13 5 13 5 

t v/ t , t / 

2 4 6 2 4^6 

Lema 3.2 Dada uma árvore, H = {V, T), de raiz v, e um vértice w & V , existe um único passeio 
que parte de v e termina em w, a saber, {v,r~^{w), . . .r~'^{w),r~^{w),w). Ademais, este passeio 
é um caminho. 

3.2 Relagoes de Ordem 

Uma ordem num conjunto S é uma relagáo, <, tal que Va, h,c & S, temos que 

1. a < a. 

2. a <h A h < c ^ a < c. 

i.e., uma relagáo reflexiva e transitiva. Alternativamente denotaremos a < h por h > a. 
Sendo = a relagáo identidade e ^ a negagáo da relagáo de ordem, dizemos que a ordem é 

• total, sse a "^ h ^ h < a. 

• parcial, sse {a < h A h < a) ^ a = h. 

• boa, se é parcial e total. 



Lema 3.3 Dado um grafo, G = (A^, F), a fungáo de descendéncia, T, define uma ordem em seus 
vértices, a ordem natural, <, definida por por j > i -^ j & f (?). Note, porém, que a ordem natural 
náo é, em geral, nem parcial nem total. 

Exemplo 3: 
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Nos grafos seguintes, de matrizes de adjacéncia Bi, B^, B^ e B^, 
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temos ordens naturais 
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• em Gi. Nem parcial, pois 1 < 2 < 1, nem total, pois 1 ^ 3 ^ 1. 

• em G^: Total, mas náo parcial, pois 1 < 4 < 1. 

• em G^: Parcial, mas náo total, pois 2 ^ 3 ^ 2. 

• em G4: Boa. 

Uma equivaléncia num conjunto S é uma relagáo, ~, tal que Va, h,c E S: 

1. a ^ a. 

2. a~6 A h ^ c^ a ^ c. 

3. a r^ h ^ h r^ a. 

i.e., uma relagáo reflexiva, transitiva e simétrica. 

A classe de equivaléncia de um elemento qualquer, a G S*, é o sub-conjunto de S: [a] = 
{x E S \ X r^ a}. 

Uma partigáo de um conjunto S é uma colecáo P de sub-conjuntos náo vazios de S tal que: 

1. yx,Y eP, X ^ r ^ X n r = 0. 

2. UxepX = S. 
i.e., uma colegáo de conjuntos disjuntos que reunidos é igual a S. 

Lema 3.4 Dado S , um conjunto munido de uma ordem, Va, hES,ar^h^a<h Ah < a, define 
uma relagáo de equivaléncia. Esta é a equivaléncia induzida pela ordem. Dado S um conjunto 
munido de uma equivaléncia, conjunto das classes de equivaléncia em S é uma partigáo de S . 

As componentes fortemente conexas, CFC, de um grafo G sáo as classes de equivaléncia 
induzida pela ordem natural nos vértices de G. G é dito fortemente conexo sse possui uma única 
CFC. 

Exemplo 4: 

As CFC dos grafos do Exemplo 3 sáo: 
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• emG'i: V^ = {{1, 2}, {3}, {4}}. 

. emG^: V = {{1,2}, {3,4}}. 

. emGseG^: V = {{1} , {2} , {3} , {4}} . 

O grafo reduzido, de um grafo G = (V,r), é o grafo G = {V,r), que tem por vértices as 
CFCs de G: V = {Vi, V^, . . . V^}, e a fungáo de filiagáo, F, definida por 

Vs G f(V;) ^ 3í G Vr, jeVslje T{{}, {Vr ^ Vs) . 

Exemplo 5: 

As matrizes de adjacéncia dos grafos reduzidos dos grafos do Exemplo 3, com os vértices 
correspondendo as CFCs na ordem em que aparecem no exemplo anterior, sáo Bi, B^, -B3 e B^. 
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Dado um conjunto S e < uma ordem em S, a ordem reduzida é a relagao de ordem no 
conjunto das classes (da equivaléncia induzida pela ordem <), S", definida por 

VX, F G ^ , X <Y ^3xeX,^3yeY\x<y . 

Lema 3.5 ^ ordem reduzida (nas CFCs) de um grafo G = {V,r), é a ordem natural do grafo 
reduzido G = {V , V), i.e., se V = {Vi, . . . Vk}, 

Vr<Vs ^ VsG t{Vr) 

^ 3ie Vr, j eVs\j e T{i) 
^ ^ieVr, j eVs\j >i . 

Ademais, 

1. A ordem natural de G é parcial sse, a menos de loops, G é acíclico. 

2. Crafos reduzidos sáo acíclicos e ordens reduzidas sao parciais. 



Teorema 3.1 (Hoffman) ; Um grafo G = {V,r) é fortemente conexo sse dado qualquer sub- 
conjunto próprio dos vértices W G V, W ^ V , houver uma aresta de um vértice em W para um 
vértice fora de W . 
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Demonstracáo: 

Se houver W um subconjunto próprio de V tal que r(W^) = 0, entáo náo há caminho de 
nenhum vértice w & W para nenhum vértice v & V, e G náo é fortemente conexo. Presupondo 
a condicáo ^W ^ V \ r(W) = 0, provemos que existe um caminho do vértice w ao vértice v, 
Vw,f G V. Tomemos inicialmente Wq = {w}. A condigáo garante a existéncia de um caminho de 
comprimento (número de arestas) 1 de w a algum vértice xi j^ w. Se xi = f a prova está completa. 
Caso contrário tomemos Wi = {^,3:1}. A condigáo garante a existéncia de um caminho c^ de 
comprimento jc^l < 2 de w para algum vértice x^ 7^ w,Xi. Repetindo o argumento até obtermos 
Xk = V, k < n, concluímos a demonstragáo, QED. 

Dado um conjunto S, munido de uma ordem <, a ordem de S, uma boa ordem, <=, em S é 
dita uma ordem coerente (com <) sse: 

1. ya,b e S, a <b A b ^ a^ a <= b . 

2. Va, b,c E S, a <= b <= c A a^^c^ar^br^c. 

O primeiro critério de coeréncia determina que o reordenamento se subordine á ordem parcial 
do grafo reduzido; O segundo critério determina que se discriminem (náo se misturem) vértices 
em CFCs incomparáveis porém distintas. 

Uma boa ordem, ou reordenamento, no conjunto A^ = {1,2, . . .n}, qi <= q^ <= . . . <= qn, 
corresponde a um vetor de permutagáo q = [qi, q^, ■ ■ ■ qn] = [o'(l)) ^"(2), . . . o'(n)]. Um reordena- 
mento coerente dos vértices de um grafo G é um reordenamento coerente com a ordem natural 
do grafo. Um reordenamento coerente num grafo acíchco é também dito um (re)ordenamento 
topológico (dos vértices) deste grafo. 

Exemplo 6: 

Listamos a seguir alguns reordenamentos nos grafos do Exemplo 3, indicando se satisfazem aos 
dois critérios de coeréncia. No grafo Gi: q = [4; 1, 2; 3] (1, 2), g = [1, 3, 2, 4] (1), g = [3, 4, 1, 2] (2), 
q = [1,3,4,2] (). No grafo G^: q = [3,4,1,2] (1,2), q = [1,3,2,4] (1). No grafo ^3: q = 
[1, 3, 2, 4] (1, 2), g = [1, 2, 3, 4] (1, 2), e estes sáo os únicos reordenamentos coerentes. No grafo G^: 
O único reordenamento coerente é q = [1, 2, 3, 4]. 

Para reordenar coerentemente os vértices de um grafo precisamos pois determinar o grafo 
reduzido, e no grafo reduzido uma ordem topológica. Esta tarefa pode ser levada a cabo com o 
algoritmo de Tarjan, a ser visto a seguir. 



3.3 Busca em Profundidade 

A busca em profundidade é um procedimento que nos permite visitar, a partir de um determinado 
vértice v de G = [N, F), todos os seus descendentes. Cada vértice será marcado como "já visitado" 
ou "náo visitado" , sendo "náo visitado" o estado inicial de todos os vértices. 
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Na busca em profundidade, partindo de v, seguiremos um caminho, sempre por vértices náo 
visitados, o mais longe possível. Ao passar por um vértice, marcá-lo-emos como "já visitado". 
Náo sendo mais possível prosseguir de um dado vértice, isto é, quando estivermos num vértice 
sem filhos náo visitados, retornaremos ao vértice de onde este foi atingido e prosseguimos na 
busca em profundidade. Quando tivermos voltado ao vértice inicial, v, e já tivermos visitado 
todos os seus filhos, teremos terminado a busca em profundidade. Os vértices visitados e as 
arestas utihzadas para atingí-los formam uma árvore de raiz v, que é a árvore da busca. 

Iniciando novas buscas em profundidade, nas quais consideramos "náo visitados" apenas os 
vértices que náo pertencem a nenhuma árvore previamente formada, teremos uma fioresta que 
cobre o grafo G, a fioresta de busca. 

Uma maneira de rotularmos (ou reordenarmos) os vértices de G durante a busca em profun- 
didade é pela ordem de visitagáo, Ov{), onde Ov{i) = k significa que o vértice i foi o /c-ésimo 
vértice a ser visitado na formagáo da fioresta de busca. Uma maneira ahernativa de rotular (re- 
ordenar) os vértices de uma fioresta de busca é a ordem de retorno: Or(), é a ordem em que 
verificamos já termos visitado todos os filhos de um vértice e retornamos ao seu pai na árvore (ou 
terminamos a árvore se se tratar de uma raiz). 

Estando os vértices de um grafo bem ordenados por algum critério, por exemplo pelos seus 
índices, a fioresta de busca em profundidade canónica é aquela em que tomamos os vértices, tanto 
para raízes de novas árvores quando para visitagáo dentro de uma busca, na ordem estabelecida 
por este critério. No exemplo 7, a primeira é a fioresta canonica. 

Exemplo 7: 

Um grafo G, e várias das possíveis fiorestas de busca que o cobrem, estáo dadas na figura 
seguinte. Apresentamos estas fiorestas com os vértices numerados primeiro pela ordem de visitagáo, 
depois pela ordem de retorno. Indicamos com um circunfiexo as raizes da busca. 
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Descrevemos agora o algoritmo de Tarjan para determinacao das componentes fortemente 
conexas de um grafo G. 

1. Considerando a boa ordem dos índices, construa a floresta de busca canónica em G, marcando 
seus vértices pela ordem de retorno; 

2. Considere G^^ com os vértices reordenados, isto é rotulados, na ordem inversa da ordem de 
retorno estabelecida no passo 1; 

3. Considerando a boa ordem estabelecida no passo 2, construa a floresta de busca canónica 
em G~^. 

Teorema 3.2 (Tarjan) ; Cada árvore em G^^ construída no passo 3 do algoritmo de Tarjan 
cobre os vértices de exatamente uma componente fortemente conexa de G. Mais ainda, a ordem 
de visitagáo (bem como a ordem de retorno), na floresta canónica do passo 3 do Algoritmo de 
Tarjan, é um reordenamento coerente dos vértices de G. 

Demonstracáo: 

Se v,w & V estáo numa mesma componente de G entáo certamente pertencem a uma mesma 
árvore em G~^ (bem como em G). Se x é um vértice de G, denotaremos o número que marca o 
vértice x, pela ordem de retorno em G, por Or{x). Se w é um vértice de uma árvore de raiz v, 
em G~^, entáo: 

1. v descende de w, em G, pois w descende de v em G~^. 

2. w descende de v, em G. 

Para justiflcar a segunda aflrmagáo notemos que, por construcáo Or{v) > Or{w). Isto signiflca 
que ou w foi visitado durante a busca em profundidade a partir de v, ou w foi visitado antes de 
V. A primeira hipótese implica em (2) enquanto a segunda hipótese é impossível pois, como por 
(1) V descende de w, jamais poderíamos ter deixado v "nao visitado" após concluir uma busca em 
profundidade que visitasse w. Q.E.D. 

Exemplo 8: 

A flgura seguinte apresenta os passos na determinagáo pelo algoritmo de Tarjan, das CFCs do 
exemplo 7, e nos dá um reordenamento coerente de seus vértices. A floresta de BEP em G está 
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rotulada pela ordem de retorno, e a floresta em G~^ pela ordem de visitacáo. O reordenamento 
coerente dos vértices do grafo original obtido neste exemplo é g = [3, 6, 7 | 1 | 4, 5, 8 | 2]. 
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3.4 Grafos Simétricos e Casamentos 

Um grafo G = {N, T) é dito simétrico sse Vi, j & N, j E T{i) =^ i E T{j). É usual representarmos 
um grafo simétrico substituindo cada par de arestas, {i,j) e {j,i), por uma aresta sem orientagáo 
ou lado, {i,j}. Podemos deflnir um grafo simétrico através de seus vértices e dos seus lados, 
G = {N, E), onde cada lado é um conjunto de dois vértices e {i,j} E E -^ i E T{j). 

Um m-casamento, num grafo simétrico, é um conjunto de m lados onde sáo distintos todos 
os vértices; M = {{ui,U2},{u3,U4}, . . . {u2m-i,U2m}}- Os vértices em M dizem-se casados, os 
de mesmo lado dizem-se companheiros, e os vértices fora do casamento dizem-se solteiros. Um 
m-casamento é máximo em G se nao houver em G um m + 1 casamento e é perfeito se náo deixar 
nenhum vértice solteiro. 

Dado um grafo simétrico G = {V, E) e nele um m-casamento, M, um caminho M-alternado 
é um caminho G = {{wo,wi}, {wi, w^}, . . . {wk^i,Wk}) que tem lados, alternadamente, dentro e 
fora do casamento. Um caminho M-ahernado diz-se um caminho de aumento se comega e 
termina em vértices solteiros. 

Teorema 3.3 (Berge) .■ Dado um grafo simétrico, G = {V, E) e nele um m-casamento, M = 
{{ui,U2}, {«3,^4}, . . . {u2m,U2m}}, estc casamcnto é máximo se náo houver caminho de aumento. 

Demonstracáo: 

Suponha que existe um caminho de aumento 
G = {{wq, wi}, {wi, W2}, . . . {W2fc-i, W2k+i}) entáo 

M' = M®G = M-{{wi,W2},{w^,Wi},...{w2k-l,UJ2k}} 
+ {{Wq, Wi}, {W2, W3}, . . . {W2k, W2k+l}} 
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é um m + 1-casamento. 

Se, por outro lado, M náo é máximo, seja M'um m+1-casamento e considere o grafo auxiliar 
H de lados E^ = M' © M, sendo Vh = Vg- Cada vértice de ií é isolado, ou pertence a algum 
lado de En e no máximo a dois lados, um de M e outro de M'. Cada componente de -fí é pois, 
ou um vértice isolado, ou um ciclo de lados alternadamente em M e M', ou um caminho de lados 
alternadamente em cada um dos casamentos. Como H tem mais lados de M'que de M, há ao 
menos uma componente de tipo caminho que comega e termina com lados de M'. Este caminho 
é M-alternado e seus vértices extremos sáo,por construgáo, solteiros em M. Temos assim um 
caminho de aumento para M. Q.E.D. 

Exemplo 9: 

Apresentamos agora: Na primeira hnha: M: um m-casamento num grafo, C: um caminho de 
aumento, e M': o m + 1-casamento M' = M ®C . Na segunda hnha: M: um m-casamento, M': 
um m + 1-casamento, e H: o grafo H = M' © M. 



4 ••• 5 ••• 6 4 5 ••• 6 4 ... 5 _ 6 

7 ••• 8 ••• 9 7 8 ••• 9 7 ... 8 - 9 
1 ••• 2 ••• 3 1-2 ••• 3 1 ••• 2 3 
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Um grafo C = {V, T) é bipartido se houver uma bipartigáo de seus vértices tal que todos os 
lados tenham um vértice em cada pedago da bipartigáo, i.e. se for possível encontrar conjuntos 

X,Y\V = XUY A Xnr = A 'ÍeeE,e = {x,y} xeX^yeY. 

Teorema 3.4 (Hall) .■ Dado C = {V., T) um grafo simétrico e bipartido, com bipartigáo V = 
X UY , existe em C um casamento que casa todos os vértices de X sse vale a condigáo de Hall: 
\/S C X, ^T{S) > i^S . isto é, qualquer subconjunto de X tem, coletivamente, ao menos tantos 
filhos quanto elementos. 

Demonstragáo 

Se existe um casamento M que casa todos os vértices em X, entáo os lados do casamento 
garantem ao menos a igualdade na condigáo de Hall. Por outro lado, se houver um casamento 
máximo, M, que deixe algum x E X solteiro, exibiremos um S* C X que viola a condigáo de HaU. 
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Seja Z o conjunto dos vértices conexos a x por caminhos M-alternados. Pelo teorema de Berge 
sabemos que náo há solteiro em Z, pois caso contrário teríamos um caminho de aumento e M náo 
seria máximo. Assim, um caminho M-ahernado maximal, isto é, que náo pode ser continuado, 
que comega no soheiro x E S (^ X, deve necessariamente terminar, com um número par de lados, 
num casado x' G 5 C X. 

Considerando, pois, S = {Z n X) + x e T = Z HY, temos que ^S = #T + 1 e T{S) = T, 
donde ^r{S) = #T = ^S — 1 < i^S, mostrando que S viola a condicáo de HalL Q.E.D. 



3.5 O Algoritmo Húngaro 

Dado um grafo simétrico G = {V,E), bipartido em conjuntos de mesma cardinahdade X eY, o 
algoritmo húngaro fornece um casamento perfeito ou encontra um conjunto S G X que viola a 
condicáo de HaU. 

Seja M um casamento que deixa x G X solteiro. Uma árvore de raiz x, H = {Vh, Eh), subgrafo 
de G, é M-ahernada se for uma árvore onde qualquer caminho partindo de x, C = {x, Wi, w^, . . .) 
for M-alternado. Observemos que se uma dada árvore M-alternada de raiz x tiver uma folha 
soheira, entáo o caminho C que vai da raiz a esta folha é um caminho de aumento e M' = M ®C 
é um novo casamento onde, além de todos os vértices casados em M, x também é casado. 

Se, todavia, encontrarmos uma árvore M-alternada de raiz solteira, onde todas as folhas sáo 
casadas e que seja máxima, isto é, tal que náo exista nenhum vértice em Vg — Vh, adjacente a 
uma folha de H teremos encontrado um S <Z X que viola a condigáo de Hall, S = Vh ^ X e 
T = VHr\Y = V{S). 

O algoritmo Húngaro faz o seguinte: dado um grafo bipartido e simétrico com um casamento 
que deixa x E X soheiro, gerar a árvore canonica de busca em profundidade M-alternada, até 
encontrar um caminho de aumento, isto é, uma folha solteira. Caso isto náo seja possível con- 
statarmos a violacáo da condigáo de HaU. 

Exemplo 10: 

No exemplo 10 temos dois exemplos de grafo e casamento, nos quais está assinalado um 
vértice w soheiro e sem mais nenhum pretendente disponível. Temos as respectivas árvores M- 
alternadas, marcadas pela ordem de visitagáo, os caminhos de aumento, C, e os novos casamentos, 

M' = M (B C, que casam o vértice w. 

O algoritmo Húngaro, bem como qualquer algoritmo casamenteiro conhecido, náo é hnear; 
Portanto vale a pena procurar heurísticas mais eficientes para a procura de um casamento perfeito. 
Dado um grafo e um casamento, os pretendentes de um vértice soheiro sáo os soheiros em seu 
conjunto de adjacéncia, e o peso de um solteiro, p{x), é seu numero de pretendentes. A idéia 
subjacente em todas estas heurísticas é a de casar primeiro os vértices mais exigentes, i.e. de 
menor peso, de modo a preservar o máximo de pretendentes para os vértices ainda solteiros. 
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Na Heurística do Mínimo Simples, HMS, dado um m-casamento imperfeito, procuraremos 
obter um m+1-casamento adicionando ao casamento um lado que una vértices ainda solteiros. 
Escolliemos este lado, de modo a casar um vértice de mínimo peso, entre os ainda solteiros, com 
um de mínimo peso dentre seus pretendentes. Sempre que houver um vértice isolado, i.e. sem 
pretendentes ou de peso zero, procuraremos casá-lo pelo algoritmo Húngaro. Na Heurística do 
Mínimo Par, HMP, escolhemos um lado a ser acrescentado ao casamento onde um dos vértices 
é de peso mínimo, e o outro tenha peso mínimo dentre todos os pretendentes a vértices de peso 
mínimo. Como a HMP pode ser bem mais custosa podemos, por exemplo, utilizar a HMS enquanto 
o peso mínimo dos vértices ainda soheiros for grande, i.e. acima de um dado hmite, e a HMP 
caso contrário. 

Exemplo 11: 

O exemplo seguinte ilustra aaphcagáo do algoritmo húngaro com a heurística min-min. No 
exemplo indicamos o peso, ou número de pretendentes, de cada vértice. Indicamos ainda, a cada 
passo, o vértice de nínimo peso a ser casado, bem como, em itálico, seu pretendente de mínimo 
peso. 
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A árvore M-alternada com raiz em x^, que ficou isolado, nos fornece o caminho de aumento 



c = X2 - yi - xs 



1/3 



de modo que podemos recomegar a heurrística em M (B c: 
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Capítulo 4 

ELIMINAgÁO ASSIMÉTRICA 



Esparsidade na Fatoragáo LU 

A maioria dos sistemas lineares encontrados na prática, especialmente os sistemas grandes, sáo 
esparsos, isto é, apenas uma pequena parte dos coeficientes do sistema náo sáo nulos. Ao tratar 
um sistema esparso, devemos tentar manter a esparsidade do sistema ao longo das tranformagóes 
necessárias á sua solugáo. Tal medida visa, primordialmente, minimizar o número de operagóes 
aritméticas a serem realizadas, na própria transformagáo e na solugáo do sistema. Economia de 
memória também é um fator importante. 



4.1 Preenchimento Local 

No método de Gauss as tranformagóes ''A -^ '^^^A podem diminuir o número de elementos nulos. 
Se tal ocorrer, dizemos que houve preenchimento de algumas posigoes. Queremos escolher o 
elemento pivó em '^A de modo a minimizar este preenchimento. 

Teorema 4.1 (Tewarson) Seja 



k M 



A[ 



k An 



"A 



k /ifc+l 



fc+1 



k An 



'A 







M^+i 



fc+l 



k An 



'A 



seja '^B, {n — k) X [n — k), a matriz Booleana associada a suhmatriz das n — k últimas linhas e 
colunas de '^A, i.e. 



'B = BiAtXi 



k+l:n\ 
+l:ra/' 



OU 



'Bl = 1 ^ '^AlXl ^ 
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e seja 

kQ^ kQ^k^y kQ _ 

onde operador complemento aplicado á matriz B, B, troca O's por l's e vice-versa. 
A escolha do pivó ^A^l ^ implica no prenchimento de exatamente ^G\ posigóes. 

Demonstragáo: 

Seja ^A a matriz obtida de ^A por permutagáo da /c + 1-ésima linha e coluna com, respectiva- 
mente, a A; + i-ésima linha e a /c + j-ésima coluna. 

Os novos elementos de '''^^A seráo, para p,q = 2 . . .{n — k), 

k+l Ak+q _ 
^k+p ~ 

k Xk+1 k+l ]\/fk+l k 7k+q 

— ^k+p ^^^k+p ^k+l 

k ^k+q k Ak+Í k ^k+q i k 'Ák+l 

— ^fc+p ^k+p ^fc+l / ^fc+l 

k /ik+s k Ak+j k /\k+s i k /\k+j 

~ ^fc+r ^fe+r ^fe+i / ^fe+í 

onde 

r = p se p ^ i, e r = 1 se p = i ; e s = gseg7^j, es = lseg = j. 

Haverá preenchimento de uma posigáo em ^^^A sempre que, 
para r, s = 1 . . . {n — k) , r ^ i e s ^ j , 

''B'^ = A ''BÍ = 1 A ''B¡ = 1 

e portanto o total de preenchimentos, correspondente ao pivo {k -\-i, k + j), é: 

n—k n—k 

J2 E ''^r ^Bi ''Bt = 

r=l, r^i s=l, S/íj 
n—k 

= Y: 'Bi{{'^B)'):'^B¡ 

r,s=l 

= { ''B{ ''B)' ''B)¡ = G{ 

Na penúltima passagem usamos que os termos r = i ou s = j sáo todos nulos, 
pois BlBl = 0. QED. 

Observagáo 4.1 Uma aproximagáo para a matriz '^G é a matriz de Markowitz: 

''FÍ = {''B1 ''B)Í , 

onde denotamos por 1 a matriz quadrada de l's. A aproximagáo F é exatamente o número de 
multiphdores náo nulos vezes o número de elementos náo nulos, fora o pivo, na hnha pivó. Esta 
aproximagáo é bastante boa se ''B é muito esparsa. 
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Exemplo 1: 

Dada a matriz ^A, cujos elementos náo nulos estáo indicados na matriz boleana associada 
B{ ^A), indique uma escolha de pivós que minimize preeenchimentos locais. Assuma que ao longo 
do processo de triangularizagáo náo há cancelamentos, i.e., que um elemento náo nulo uma vez 
preenchido, náo voha a se anular. 



Tomando 



^B = B{ M) = °5 



1110 
110 
10 10 
11 
1110 1 



^G 



3 15 3 
10 3 2 

2 3 3 12 

3 6 5 11 
16 3 6 3 



Assim, arg min 



{i,j)\''B¡=l 



Lt,- 



{(1,1), (2, 3)}. 



Escolhendo {i,i) = (1, 1), i.e. A\ como pivó, temos o preenchimento de G\ = posigoes em 



U, 



^B = B{ ^A+ ^M) 
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11 
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2 



Assim, ar5(min(.^^.-)|i^,^-^ ^Gl ={(1,2)}. 

Escolhendo {i,j) = (1,2), i.e. ^A\ como pivó, temos o preenchimento de ^G\ = posigóes em 



'A, 



''B 



Bi^A 



'M) 



1110 
110 
110 
11 
i i 1 1 



'G 



1 

2 
1 



1 
1 
2 



Assim, argmm^. .^^2^1=1 



Escolhendo {i,j) 



'G¡ ={(1,1), (1,2), (2, 2), (2, 3), (3,1), (3, 3)}. 
1, 1), i.e. ^Ag como pivo, temos o preenchimento de ^Gj^ = 1 posigóes em 



'A, 



''B = B{ ^A+ ^M) 



1 


1 


1 
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1 














1 


1 














1 


1 


1 


1 


1 


1 


1 



'B 



1 1 
1 1 



Como ^B = 1, é obvio que nao haverá mais preenchimento, quaisquer que que sejam os pivós 



doravante selecionados. Tomando, por exemplo, o pivo Ag, temos finalmente temos a fatoragao 
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A = PAQ = LU, com 



B{Á) 
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B{M + U) 



1110 
110 
110 

11111 

11 



sendo os vetores de indices permutados *^p = [1,2,3,5,4] e "^g = [1,3,2,4,5]. Neste processo de 
triangulaziragáo, tivemos assim o preenchimento de apenas 1 posigáo, em ^A\. // 

Observagáo 4.2 Para náo fazer a selegáo de pivós apenas pelos critérios para minimizagáo de 
preenchimento local, desprezando assim completamente os apectos de estabihdade numérica, pode- 
mos "vetar" escolhas de pivós muito pequenos, i.é., que imphquem no uso de muhiphcadores muito 
grandes, \M¡\ > multmax ^ 1. Se nosso sistema for bem equihbrado, podemos usar como critério 
de veto o tamanho do proprio pivó, \Ul\ < pivomin <^ 1. 



4.2 Pré-Posicionamento de Pivós 



Minimizar o preenchimento local é um processo custoso, tanto no número de operagóes envolvidas 
no cálculo de ''G, como por termos de acessar toda a matriz ''A, o que é um grande inconveniente 
para o armazenamento eficiente da matriz. Ademais, minimizar o preenchimento local, i.e. aquele 
causado por cada pivo individualmente, é uma estratégia gulosa que pode ter um fraco desem- 
penho global, i.e. causar muito mais preenchimento que a seqüéncia ótima de n pivos. Visando 
superar as deficiéncias dos métodos locais estudaremos a heurística P3, ou Procedimento de 
Pré-determinagáo de Pivós. 

A heurística P3 procura uma permutagáo de hnhas e colunas, B = PAQ, que reduza o preenchi- 
mento na fatoragáo B = LU. O procedimento que implementa a heurística P3 posiciona as colunas 
de A em B e numa matriz temporária, S, e permuta hnhas de todas as colunas, tentando produzir 
uma B que seja quase triangular inferior. As colunas de B que náo tem ENNs acima da diagonal 
denominam-se colunas triangulares. Colunas triangulares sáo sempre posicionadas com um 
ENN na diagonal, que será usado como pivo da coluna na fatoragáo B = LU. As colunas de 
B que tem ENNs acima da diagonal denominam-se espinhos. Na fatoragáo LU de B, somente 
poderá ocorrer preenchimento dentro dos espinhos. Ademais, náo haverá preenchimento acima do 
mais aho (menor índice de hnha) ENN num espinho. Portanto, queremos minimizar o número de 
espinhos e, como objetivo secundário, minimizar a ahura dos espinhos acima da diagonaL 

A heurística P3 procede como segue: 

1. Compute o número de ENNs em cada hnha e coluna de A, ou os pesos de hnha e coluna, 
respectivamente, pr{i) e pc{j). 
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2. Compute h = argminipr{i) e p = pr{h). 

3. Compute a p-altura de cada coluna, Qp{j), 

o numero de ENNs na coluna j em linhas de peso p. 

4. Compute t = argmaxj Qp{j). 

(a) Se p = 1, seja h E {i \ pr{i) = 1 A A{i,t) ^ 0}; posicione t como a primeira coluna 
de A, h como a primeira linha, e aphque P3 recursivamente a y4(2:m,2:n). A coluna 
excluída torna-se a úhima coluna de B. 

(b) Se p > 1, posicione a coluna t como a úhima coluna de A, e aphque P3 recursivamente 
a A(:,l:n-1). A coluna excluída torna-se a primeira coluna de S. 

(c) Se p = 0, posicione h como a primeira hnha, reposicione a primeira coluna de S* (úhima 
a ser excluída de Á) como a úhima coluna de B, e aphque P3 recursivamente a A(2:m,:). 

Em P3, o caso 

• p = 1 corresponde ao posicionamento de uma coluna triangular em B. 

• p > 1 corresponde a impossibihdade de posicionar uma coluna triangular. Assim (tempo- 
rariamente) ehminamos uma coluna de A e prosseguimos com P3. O critério de selegáo para 
a coluna a ehminar visa produzir o caso p = 1 nos próximos passos de P3. Em caso de 
empate na p-ahura, poderíamos usar como desempate a (p -|- l)-ahura. 

• p = corresponde a nao haver em A uma coluna que pudéssemos posicionar em B de modo 
a continuar formando B náo singular. Entáo reintroduzimos, como próxima coluna de B, a. 
primeira coluna em S. Escolhemos a primeira em S visando minimizar a ahura do espinho 
acima da diagonal. Existe o perigo de que este espinho náo venha a prover um pivó náo nulo, 
ou que após cancelamentos, o pivo seja muito pequeno para garantir a estabihdade numérica 
da fatoracáo. Se A é esparsa e tivermos apenas alguns espinhos, a melhor solugáo é intercalar 
a fatoragáo numérica e simbóhca, o que nos dá a chance de rejeitar pivós instáveis. Métodos 
para hdar com pivós inaceitáveis na fase de fatoragáo numérica, após termos completado 
uma fase independente de fatoragáo simbóhca, sáo discutidos posteriormente. 

Exemplo 

Aphquemos o P3 á matriz booleana A. Á medida que o P3 prossegue representaremos a matriz 
particionada B, A, S , cujos ENNs seráo denotados respectivamente b, ae s. Os índices de hnha 
e coluna estáo a esquerda e acima da matriz. Os pesos, p, e as p-ahuras, Op, a direita e abaixo 
da matriz. Os pivós, ou espinhos, escolhidos a cada passo estáo em negrito. 
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Na permutagáo final obtida pela P3, dada pelos vetores de permutagáo p e q, indicamos com um 
+ as posiccoes a serem preenchidas no processo de eliminagáo. Uma notagáo mais compacta para 
o mesmo exemplo seria: 
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No exemplo seguinte reintroduzimos um espinho com na posigao pivó. Todavia, no processo 
de eliminagáo, esta posigao será preenchida antes da sua utihzagáo como pivó. 
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Observagáo 4.3 Da maneira como descrevemos o P3, poderíamos aphcá-lo também a uma matriz 
retangular. No caso de Programagáo hnear é comum ordenarmos pelo P3 todas as colunas da 
matriz de restrigóes, A, e depois tomarmos, a cada reinversáo, as colunas na base, B, conforme a 
ordem estabelecida por P3 em A [Orchard-Hays68] . 
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Exercícios 

1. Implemente a fatoracáo LU usando a heurística de Markowitz. Antes de aceitar um pivó, 
assegure-se que este tem módulo maior que pivmin = 10""', substituindo-o caso contrário. 
Use a representagáo estática por colunas e de rede da matriz. 

2. Implemente o P3 para ordenar as colunas de uma matriz retangular, conforme a observacáo 
6.3. Use a representagáo estática por coluna da matriz. 

3. Considere a possibilidade de termos um espinho problemático que, após o processo de elim- 
inagáo das colunas precedentes, apresente um na posigáo pivo. Mostre que necessariamente 
existe um espinho, á direita do problemático, cujo elemento na mesma linha do pivo do es- 
pinho problemático, neste estágio da fatoragáo da matriz, é diferente de zero. Descreva 
um procedimento para lidar com estes casos excepcionais através da permutagáo de colunas 
espinhos. Discuta a viabilidade de resolver o problema através de permutagáo de linhas. 



Capítulo 5 

FATORAgÓES SIMÉTRICAS e 
ORTOGONAIS 

Projetores e Problemas Quadráticos 
5.1 Matrizes Ortogonais 

Dizemos que uma matriz quadrada e real é ortogonal sse sua transposta é igual a sua inversa.Dada 
Q uma matriz ortogonal, suas colunas formam uma base ortonormal de 3?", como pode ser visto 
da identidade Q'Q = 1.0 quadrado da norma quadrática de um vetor v, 



n 
I l|2 _ V^/ \2 / 



permanece inalterada por uma transformacao ortogonal, pois 

\\Qv\\ = {Qv)'{Qv) = v'Q'Qv = v'Iv = v'v = \\v\\ . 



5.2 Fatoragáo QR 



Dada uma matriz real A m x n, m > n, podemos existe uma matriz ortogonal Q tal que A 

D 

, onde R é uma matriz quadrada e triangular superior. Esta decomposigáo é dita uma 



Q 







fatoracao QR, ou fatoragáo ortogonal, da matriz A. Descrevemos a seguir um método para 
fatoracáo ortogonal. 



A rotagáo de um vetor 



Xi 
X2 



G 3?^ por um ángulo 6 é dada pela transformacao linear 
rot{6)x - 



cos{9) sin(í^) 
— sin(í^) cos(í^) 



Xi 
X2 
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Notemos que a rotagao é uma transformagao ortogonal, pois 



rot{eyrot{e) 



cos(^)^ + sin(^)^ 






cos(^)^ + sin(^)^ 



1 
1 



Além disso podemos tomar o ángulo e = — arctan(a;2/a;i) de modo que a rotacao correspon- 
dente anule a segunda componente do vetor rodado (se xi = 0, tome e = 7r/2). 

Como o produto de transformagoes ortogonais continua ortogonal (prove), podemos usar uma 
seqüéncia de rotagóes para levar a matriz A á forma triangular superior, como veremos a seguir. 

A rotagáo de Givens é um operador linear cuja matriz coincide com a identidade, exceto 
num par de linhas onde imergimos uma matriz de rotagáo bidimencional: 



G{t,j,e) 



cos(^) 



sin(^) 



— sin(^) cos(í^) 



Dizemos que a aplicagao deste operador numa matriz A, GA, roda as linhas i e j de A de um 
ángulo e. 

Abaixo ilustramos uma seqüéncia de rotagóes de hnhas que leva uma matriz 5 x 3 á forma 
triangular superior. Cada par de índices, {i,j), indica que rodamos estas hnhas do ángulo apro- 
priado para zerar a posigáo na hnha i, coluna j. Supomos que, inicialmente, a matriz é densa, i.e. 
todos os seus elementos sáo diferentes de zero, e ilustramos o padráo de esparsidade da matriz nos 
estágios assinalados com um asterisco na seqüéncia de rotagóes. 

(1,5)*(1,4)(1,3)(1,2)*(2,5)(2,4)(2,3)*(3,5)(3,4)* 



X 


X 


X 




X 


X 


X 




X 


X 


X 




X 


X 


X 







X 


X 





X 


X 


X 







X 


X 







X 


X 







X 


X 







X 


X 





X 


X 


X 







X 


X 










X 










X 










X 





X 


X 


X 





X 


X 








X 

















_ 



Tomando Q como a produtória das rotagoes de Givens, temos a fatoragao A = QR, como procu- 
rada. 
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5.3 Espacos Vetoriais com Produto Interno 

Dados dois vetores x, y G 3fí", o seu produto escalar é definido como 

n 

< X \ y >= x'y = ^ Xiy\ 

i=l 

Com esta definicáo, o produto escalar é um operador que satisfaz as propriedades fundamentais 
de produto interno, a saber: 

1. < X \y >=< y \ X >, simetria. 

2. < ax + (3y \ z >= a < x \ z > +(3 < y \ z >, linearidade. 

3. < a; I a; >> , nao negatividade. 

4. < X \ X >= <í=^ a; = , definitividade. 



Através do produto interno, definimos a norma: 

II II I 1/2 

||x|| =< X \ X > ' ; 
e definimos também o ángulo entre dois vetores náo nulos: 

Q{x,y) = arccos(< x \y > /||a;|| ||í/||). 

5.4 Projetores 

Consideremos o subespaco linear gerado pelas colunas de uma matriz A, m x n, m > n: 

C{A) = {y = Ax,xeW}. 

Denominamos C{A) de imagem de A, e o complemento de C{A), N{A'), de espago nulo de A', 

N{A') = {y I A'y = 0}. 

O fator ortogonal Q = [C \ N] nos dá uma base ortonormal de 3?™ onde as n primeiras colunas 
sáo uma base ortonormal de C{A), e as m — n últimas colunas sáo uma base de N{A'), como pode 

ser visto diretamente da identidade Q'A = 

^ [ 

Definimos a projegáo de um vetor b G 3?™ no espaco das colunas de A, pelas relagóes: 

y = PciA)h ^ye C{A) A{b-y)± C{A) 
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ou, equivalentemente, 



y = Pc{A)b ^3x \y = Ax A A'{b - y) = 0. 



No que se segue suporemos que A tem posto pleno, i.e. que suas colunas sao linearmente 
independentes. Provemos que o projetor de b em C{A) é dado pela aplicagáo linear 

Pa = A{A'A)'^A'. 

Se y = A{{A'A)-^A'b), entáo obviamente y G C{A). Por outro lado, A'{b - y) = A'{I - 
A^A'A^-^A'^b = {A' - IA')b = 0. 



5.5 Quadrados Mínimos 

Dado um sistema superdeterminado, Ax = b onde a matriz A m x n tem m > n, dizemos que x* 
"resolve" o sistema no sentido dos quadrados mínimos, ou que x* é & "solugáo" de quadrados 
mínimos, sse x* minimiza a norma quadrática do resíduo, 

X* = Arq min \\Ax — 611, 

Dizemos também que y = Ax* é a melhor aproximagáo, no sentido dos quadrados mínimos de b 
em C{A). 

Como a multiplicagáo por uma matriz ortogonal deixa inalterada a norma quadrática de um 
vetor, podemos procurar a solugáo deste sistema (no sentido dos quadrados mínimos) minimizando 
a transformagáo ortogonal do resíduo usada na fatoragáo QR de A, 



\\Q'{Ax-h)\\' 



' R ' 


X — 


c 


_ 




. d _ 



||i?x-cf + ||Ox-d|p 



Da última expressao vé-se que a solugao, a aproximagao e o resíduo do problema original sao 
dados, respectivamente, por 



X 



R ^c , y = Ax* e z = Q 





d 



Como já havíamos observado, as m — n últimas colunas de Q formam uma base ortonormal de 
N{A'), logo z _L C{A), de modo que concluímos que y = PAbl 



5.6 Programagáo Quadrática 



O problema de programagao quadrática consiste em minimizar a fungao 

f{y) = {l/2)y'Wy + c'y, W = W' 



5.7. FATORAQAO DE CHOLESKY 
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sujeitos ás restrÍQoes 

9i{y) = N-y = di. 

Os gradientes de f e gi sáo dados, respectivamente, por 

Vyf = y'W + c' , eVyg. = N',. 

As condÍQoes de otimalidade de primeira ordem (condigóes de Lagrange) estabelecem que as 
restrigóes sejam obedecidas, e que o gradiente da fungáo sendo minimizada seja uma combinagáo 
linear dos gradientes das restrigóes. Assim a solugáo pode ser obtida em fungáo do multiplicador 
de Lagrange, i.e. do vetor / de coeficientes desta combinagáo linear, como 



N'y = d A y'W + c' = l'N' 



" iv' ■ 

_W N _ 




y 
_ i 


= 


' d' 

c 



ou em forma matricial. 



Este sistema de equagoes é conhecido como o sistema normaL O sistema normal tem por matriz 
de coeficientes uma matriz simétrica. Se a forma quadrática W for positiva definida, i.e.se 
Vx x'Wx > A x'Wx = <í=^ a; = 0, e as restrigóes A^ forem lineramente independentes, a 
matriz de coeficientes do sistema normal será também positiva definida. Estudaremos a seguir 
como adaptar a fatoragáo de Gauss a matrizes simétricas e positivas definidas. 



5.7 Fatoracáo de Cholesky 



Após a fatoragáo de Gauss de uma matriz simétrica, S = LU, podemos por em evidéncia os 
elementos diagonais de U obtendo S = LDL'. Se S for positiva definida assim o será D, de modo 
que podemos escrever D = D^^'^D^^'^, D^/"^ a matriz diagonal contendo a raiz dos elementos em 
D. Definindo C = LD^^'^, temos S = CC', a fatoragáo de Cholesky de S. 

Anahsemos agora algumas relagoes entre as fatoragóes de Cholesky e ortogonal (QR), e de que 
maneiras a fatoragáo ortogonal nos dá uma representagáo da inversa. Primeiramente observemos 
que se A = QR, 

A'A = {QR)'QR = R'Q'QR = R'IR = L'L 

i.e., o fator triangular da fatoragáo ortogonal é o transposto do fator de Cholesky da matriz 
simétrica A'A. 

Veremos agora que, no caso de termos A quadrada e de posto pleno, o produto pela inversa, 
y = A~^x, pode ser calculado sem o uso exph'cito do fator Q: Transpondo AR~^ = Q obtemos 
Q' = R-'A', donde 

y = A-'^x = {QR)-^x = R-^Q'x = R-^R-^A'x . 
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Exercícios 

1. Use as propriedades fundamentais do produto interno para provar: 

(a) A desigualdade de Cauchy-Scwartz: | < a; | y > | < ||a;|| ||í/||. Sugestao: Calcule 
ll^ — ay\\ para a =< a; | y >^ /||y||- 

(b) A Desigualdade Triangular: ||a; + í/|| < ||a;|| + \\y\\. 

(c) Em que caso temos igualdade na desigualdade de Cauchy-Schwartz? Relacione sua 
resposta com a definigáo de ángulo entre vetores. 

2. Use a definigáo do produto interno em 3Í" para provar a Lei do Paralelogramo: ||a; + y\\ + 

11 ii9 11 ii9 II ii9 

If ~ y\\ = ^IfII + ^\\y\\ ■ 

3. Uma matriz P é idempotente, ou um projetor náo ortogonal, sse P^ = P. Se P é idem- 
potente prove que: 

(a) R = {I — P) é idenpotente. 

(b) 3?" = C(P) + C(/2). 

(c) Todos os autovalores de P sáo ou +1. Sugestáo: Mostre que se é uma raíz do 
pohnómio característico de P, v^p(A) = det(P — XI), entáo (1 — A) = 1 é raíz de (/9^(A). 

4. Prove que VP idempotente e simétrico, P = Pc{p)- Sugestáo: Mostre que P'(/ — P) = 0. 

5. Prove que o operador de projegáo num dado sub-espago vetorial V, Py, é único e simétrico. 

6. Prove o teorema de Pitágoras: V6 G 3?™,« G V temos que ||6 — «11 = ||6 — PyfeH + 

II o / II 2 

\\Pvb — u\\ . 

7. Formule o problema de quadrados mínimos como um problema de programagáo quadrática. 

(a) Assuma dada uma base A^ de N{A'). 

(b) Calcule diretamente o resíduo, z = b — y, em fungáo de A. 

8. O trago de uma matriz A é definido por tr{A) = J2iA^. Mostre que 

(a) Se A, m X n, tem posto pleno p{A) = n, entáo tr{PA) = n. 

(b) Nas condigóes do item anterior, definindo Ra = {I — Pa), temos que tr{RA) = m — n. 

9. Método de Hauseholder: A reflexáo definida por um vetor unitário u \ u'u = 1, é a 
transformagáo hnear H = I — 2uu'. 

(a) Interprete geometricamente a opereagáo de refiexao. 

(b) Prove que H = H', H' = H-\ e H^ = I. 
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(c) Dado X T^ um vetor em i?", tome 



V = X ± \\x\ 



1 






■u = f /||ü|| e H = I — 2uu' 



Mostre que (Hx)i = \\x\\, e que todas as demais componentes de Hx se anulam. 
(d) Discuta como poderiamos usar uma série de reflexóes para obter a fatoracáo QR de 



uma matriz. 
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Capítulo 6 

ELIMINAgÁO SIMÉTRICA 



Esparsidade na Fatoragáo de Cholesky. 



6.1 Grafos de Eliminagáo 

Na eliminagáo assimétrica, estudada no capítulo 4, procuramos uma permutagáo geral, QPAQ' 
que minimizasse o preenchimento durante a fatoragáo QPAQ' = LU . No caso simétrico queremos 



preservar a simetria da matriz, e restringir-nos-emos a permutagoes simétricas, QAQ' 
LL'. 



A 



qU) 



Exemplo 1: 

Neste exemplo mostramos as posigóes preenchidas na fatoragáo de Cholesky de uma matriz 
simétrica A, bem como o preenchimento na fatoragáo de duas permutagóes simétricas da mesma 
matriz: 

\ X X X 

X ?) Q X X 

a; 6 a; 

X a; 2 

a; 4 

X 5 



1 
2 
3 

4 
5 
6 



1 X X 

X 2 X 

JL JL d 

X 



X 





5 X 
X X 6 



X 

4 





" 5 


4 


2 


X 


X 
X 


X 




X 






6 




X 






X 


X 




3 


X 








X 


X 


X 


1 



Assim como no capítulo 4 a hnguagem de teoria de grafos foi útil para descrever o processo 
de ehminagáo simétrica, usaremos agora grafos simétricos para estudar a ehminagáo simétrica. 
O primeiro destes conceito a ser definido é o de grafos de eliminagáo, que nada mais sáo que 
os grafos que tem por matriz de adjacéncia as submatrizes ^A do processo de fatoragáo (veja 
capítulo 2): Dado um grafo simétrico G = {N, E), N = {1, 2, . . . n}, e uma ordem de ehminagáo 
q = [o"(l), . . . o"(n)], onde a é uma permutagáo de [1, 2, . . .n], definimos o processo de ehminagáo 
dos vértices de G na ordem q como a seqüéncia de grafos de ehminagáo Gi = {Ni, Ei) onde, para 
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i = 1 . . .n, 



"5 



^i = {q'(0? Q'(^ + 1)) • • • Q'(^)}? Ei = E, e, para i > 1 , 
{a,b}eEi^ {a,b} e Ei^i ou {q{i-l),a}, {q{i - 1) , b} e Ei^i 

Definimos também o inverso ordem de eliminagáo, q{i) = k '^ q{k) = i, significando que i foi o 
/c-ésimo vértice de G a ser eUminado. 

O grafo preenchido é o grafo P = {N,F), onde F = U"i?j. Os lados originais e os lados 
preenchidos em F sáo, respectivamente, os lados em E e em F — E. 

Observagáo 6.1 Ao eUminar a j-ésima coluna na fatoragáo de Cholesky da matriz QAQ' = 
Af-/J = LL' preenchemos exatamente as posigóes correspondentes aos lados preenchidos em F 
quando da ehminagáo do vértice q{j)- 

Exemplo 2: 

Os grafos de eUminagáo e o grafo preenchido correspondentes á segunda ordem de ehminagáo 
do exemplo 1, g = [1, 3, 6, 2,4, 5], sáo: 



1-3 3 1-3 

|x\ /l\2-62 |x|\ 

2 6|2-6| x| |\ 5-4 2-61 

// //5 45-4 |x|/ 

5 4 5 4 5-4 

Lema 6.1 (Parter) Se f = {i,j} G F, entáo ou f é um lado original, i.é. F E E, ou f foi 
preenchido quando da eliminagáo de um vértice k \ q{k) < mm{q{i) , q{j)} . 



Lema 6.2 (do caminho) Considere a eliminagáo dos vértices de G = {N,E) na ordem q. 
Temos que {i,j} G F sse existe um caminho de i a j em G, passando apenas por vértices elimi- 
nados antes de i ou j , i.e., 3C = {i, Vi, v^, ■ ■ ■ Vp,j), em G, com q{l) . . . q{p) < mm{q{i), q{j)}. 

Demonstragáo: 

<í=: trivial. 

=^: por apUcagáo recursiva do lema de Parter. 

Dada a matriz A, G = {N,E) = {N,B{A)), a ordem de eUminagáo q, e o respectivo grafo 
preenchido, consideremos o conjunto de índices de Unha de ENNs na coluna j do fator de Cholesky, 
LJ I QAQ' = LL': 

enn{U) = {i\i> j ^ {q{i), q{j)} G F} + {j} . 
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Definimos a árvore de eliminaQao, H, por 

T^-i/-\ í Í; se enn(P) = {j}, ou ^ , . 

TJ(j) = < ' , / ■.. , caso contrario 

[ minji > j \ i G enn[U)\ 

Para náo sobrecarregar a notacáo usaremos /i() = T~^(y ) e 5f() = Th{ )■ Assim /i(j), o pai de j 
em ií, é simplesmente o primeiro ENN (náo diagonal) na coluna j de L. 

Exemplo 3: As as árvores de eUminagáo correspondentes ao exemplo 1 sáo: 

6^5 6 ^ 5 ^ 4 /2 

^ 1^2^3 \4-^l 

Teorema 6.1 (da árvore de eliminagáo) Dado i > j , 

i G enn{L^) =^ j E g(i) , 

i.e., qualquer mdice de linha abaixo da diagonal na coluna j de L é um ascendente de de j na 
árvore de eliminagáo. 

Demonstragáo: 



J 

X ... k 



X 



n 



Se i = h{j) o resultado é trivial. Caso contrário (vide figura 1) seja k = h{j). Mas L¿ ^ 
OAL^ ^ ^ Lf ^ 0, pois {qÍj),qii)},{(lU),QÍk)} G Gj => {q{k),q{{)} E G.+i. Agora, ou 
i = h{k), ou apUcamos o argumento recursivamente, reconstruindo o ramo de H, {i,l, . . .k,j), 
i> l > ... > k > j. QED. 

Pela demonstragáo do teorema vemos que a árvore de eliminagáo retrata as dependéncias entre 
as colunas para o processo de fatoragáo numérica da matriz. Mais exatamente, podemos eUminar 
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a coluna j de A (i.e. calcular todos os multiplicadores na coluna j, denotados por M^ no capítulo 
2, e atualizar os elementos afetados por estes multiplicadores) sse já tivermos eliminado todos os 
descendentes de j na árvore de eliminagáo. 

Se pudermos realizar processamento paralelo (veja capítulo 8), podemos eliminar simultanea- 
mente todas as colunas em um mesmo nível da árvore de eliminacáo, comegando pelas folhas, e 
terminando por eliminar a raiz. 

Exemplo 4: 

Consideremos a eliminagáo de uma matriz com o mesmo padráo de esparsidade da última 
permutagáo do exemplo 1. Sua árvore de eliminagáo é a última apresentada no exemplo 3. Esta 
árvore tem 3 níveis que, das folhas para a raiz sáo: {1, 3, 2}, {4, 5}, e {6}. Assim, podemos fatorar 
uma matriz com este padráo de esparsidade em apenas 3 etapas, como ilustrado no exemplo 
numérico seguinte: 



1 7 

2 8 

3 6 9 

7 53 2 

8 6 49 23 

9 2 23 39 



1 7 

2 8 

3 6 9 

7 4 2 

4 2 5 5 

<? 2 5 12 



1 



7 



3 6 9 

7 4 2 

4 2 5 5 

3^16 



O próximo teorema mostra uma forma computacionalmente mais eficiente de obter o grafo 
preenchido, P = {N,F), e a árvore de ehminagáo, H, dado o padráo de esparsidade da matriz 
original, G = {N,E), e a ordem de ehminagáo, q. Nesta versao simphficada dos grafos de ehm- 
inagáo, G*, ao ehminarmos o vértice q{j), preenchemos apenas os lados incidentes ao seu vizinho 
mais próximo de ser ehminado. 

Teorema 6.2 (da fatoragáo simbólica) Definimos agora a versáo simplificada do processo de 
eliminagáo: 






{N,,E*), 
E¡ = E, 
h*{j) = mm{t > j I {q{j),q{t)} G E*} 
E*+i = {{a,h}eE*\q{a),q{h)>j} 

U {{q{h{j)),v}, v\q{v)>jA{q{j),v}eE*} 
F* = \J1E* 



Afirmamos que o grafo preenchido e a árvore de eliminagáo obtidos no processo simplificado co- 
incidem com a definigáo anterior, i.e., F* = F e h* = h. 



Exemplo 5: Os grafos de ehminagao simphficados e o grafo preenchido referentes a segunda 
ordem de ehminagáo no exemplo l, q = [1, 3, 6, 2, 4, 5], sáo: 
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1-3 3 1-3 

|x \ / \ \ 2 - 6 2 |x|\ 

2 6|2 6 \ /||\ 5-42-61 

/ / //545^ I X I / 

5 4 5 4 5-4 

O lema seguinte demonstra o teorema da fatoragáo simbólica: 

Lema 6.3 

enn{U) = VJkeg(j)enn{Ü') U enn{A^) - J + {j} , J = {1, 2, . . . j} . 

Demonstracáo: 

D: 

1 

k — X 

X ... j 



n 

Consideremos k G g{j) (vide figura 2). Se i > j E enn{L^), entáo se Lj já náo é um ENN, será 
preenchido ao eliminarmos L^. 

C: 

1 



X ... h{l) 

X ... h\l) 



h^ : h^ : ■ ■ ■ : 

n 



64 CAPITULO 6. ELIMINAC^ÁO SIMETRICA 

Seja / < j' I j' G enn{Ü), i.e., consideremos uma coluna e cuja eliminacáo poderia causar preenchi- 
mentos na coluna j (vide figura 3). Pelo teorema da árvore de eliminagáo, j é um ascendente de 
/, i.e., 3p < n I j = y^^{l). Mas pela primeira parte da prova, 

enn{Ü) -JC enn{Ü'^^^) - J C . . . C enn{L''^^^^) C enn{V) - {j} , 

de modo que qualquer vértice que poderia causar, durante a eliminagáo da coluna /, uma preenchi- 
mento na coluna j, deve necessariamente aparecer em h^{l) G g{j). QED. 

6.2 Grafos Cordais 

Em um dado grafo simétrico, G = {N,E), definimos os seguintes termos: G = {N,E) é cordal 
sse para qualquer ciclo C = {vi,V2, ■ ■ ■ Vp, vi), p > 3, existe uma corda, i.e., um lado e E E hgando 
dois vértices náo consecutivos em C. Um vértice é simplicial sse sua ehminagáo nao causa 
preenchimento. Uma ordem de ehminagáo, q, é perfeita sse a ehminagáo de nenhum dos vértices, 
na ordem q, causa preenchimento. Um subconjunto dos vértices, S G N, é um separador entre 
os vértices a e b, sse a remogáo de S deixa a e b em componentes conexas distintas. Dizemos que 
S G N é um separador entre A G N e B G N sse, \/a G A, b G B, S separa a de 6. Um conjunto, 
C, satisfazendo uma propriedade, P, é minimal (em relagáo a P) sse nenhum subconjunto próprio 
de C satisfaz P. Analogamente, um conjunto é maximal, em relagáo a P, sse nenhum conjunto 
que o contenha propriamente satisfaz P. Um clique é um conjunto de vértices C G N onde todos 
os vértices sao adjacentes, i.e., Vi, j G C, {i,i} G E. indexChque 

Exemplo 6: 

No lo grafo do exemplo 2 vemos que 5 e 6 sáo vértices simphciais, q = [5,4,2,6,3, 1] é uma 
ordem de ehminagáo perfeita. No úhimo grafo do exemplo 2 vemos que S = {2,3,6,4} é um 
separador entre os vértices a = 1 e 6 = 5 (náo minimal), S' = S — {3} é um separador minimal, e 
C = {1, 2, 3, 6} é um chque. 

Teorema 6.3 (da caracterizagáo de grafos cordais) Dado G = {N,E), as trés propriedades 
seguintes sáo equivalentes: 

1. Existe em G uma ordem de eliminagáo perfeita. 

2. G é cordal. 

3. Se S é um separador minimal entre a,b G N, entáo S é um clique. 

Demonstragáo: 
1 ^2: 
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Seja q = [cr{l), cr{2), . . . cr{n)] umaordem deeliminagáo perfeitaem G, e sejaC = (fi,f2, • • .Vp,Vi] 
um ciclo em G. Consideremos Vk o primeiro vértice de C a ser eliminado, k = arg mini<j<p{g(f¿)}. 
Como q é uma ordem de eliminacáo perfeita, a corda {vk^i,Vk+i} G E. 

2^3: 

Seja S um separador minimal entre AeB,a&A, b^B,v,w&S. Como S é minimal, existe 
um caminho G = {a,Ci, . . .Cp,v) \ C\, . . .Cp^ A, pois, caso contrário, S — v continuaria separando 
a de h. Analogamente existe um caminho D conectando a &. w, com os vértices intermediários 
todos em A. Existe pois um caminho de f a w, cujos vértices intermediários estáo todos em 
A. Seja P um tal caminho de comprimento mínimo, P = {v,ai, . . .ap,w). Analogamente seja 
Q = {w, bi, . . . bq, v) um caminho de w a v através de B com comprimento mínimo. 

Concatenando P e Q obtemos um ciclo R = {v, ai, . . . a^, w, bi, . . . bq, v). Como G é cordal, R 
contém ao menos uma corda, {x,y}. Como S é um separador, náo podemos ter x & A e y E B. 
Como P tem comprimento mínimo, náo podemos ter x,y E P e, analogamente náo podemos ter 
x,y E Q. Assim, {v,w} é a única corda possível em R, e concluímos que \/v,w G S, {v,w} G G. 

Antes de 3 ^ 1, provemos 2 lemas auxihares: 

Lema 6.4 (hereditariedade) Se G satisfaz a propriedade 3, entao qualquer subgrafo de G in- 
duzido por um subconjunto de vértices também a satisfaz. 

Demonstracáo: 

Seja G = {N, E) o subgrafo de G induzido por N G N, e S minimal em G separando a de b. 
Podemos, em G, completar S, com vértices de iV — A^, de modo a obter S, um separador minimal 
em G entre a e b. Mas se G satisfaz a propriedade 3, S é um chque, e como S é um subgrafo de 
S, S também é um chque. 

Lema 6.5 (Gavril) Se G satisfaz a propriedade 3 entáo ou G é completo, ou G tem ao menos 
dois vértices simpliciais náo adjacentes. 

Demonstracáo: 

Provemos o lema por indugáo no número de vértices de G, n. Para n = 1 o lema é triviaL Se 

G = {N, E), n > 1, náo é completo, 3a, 6 G G | {a, b} ^ E. 

Seja S um separador minimal entre a e b, partindo G — S em ao menos 2 componentes conexas 
distintas, A e B, a E A, b E B. G{SUA), o subgrafo induzido pelos vértices de S U A, hereditari- 
amente satisfaz a propriedade 3, e pela hipótese de indugáo, ou G{S U A) é um chque, ou contém 
(ao menos) 2 vértices simphciais náo adjacentes. 

Se G{S U A) nao for completo, entáo um dos seus 2 vértices simphciais náo adjacentes deve 
estar em A (pois S é um chque). Se G{SUA) é completo, entáo qualquer vértice de A é simphcial 
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(pois S é um separador). Logo, existe ao menos um vértice simplicial, v E A. Analogamente, 
existe um vértice simplicial w E B, e v náo é adjacente a w (pois S separa v de w). 

3^1: 

Seja G = {N, E) satisfazendo a propriedade 3. Se G é completo, qualquer ordem q é perfeita. 
Caso contrário, existem 2 vértices simpliciais náo adjacentes. Podemos pois eliminar um deles, 
g(l), sem preenchimento, e pelo lema da hereditariedade o subgrafo resuhante continua satis- 
fazendo a propriedade 3. Podemos entáo usar o argumento recursivamente para obter uma ordem 
de ehminagáo perfeita. QED. 



6.3 Ordenacoes por Dissecgáo 

Consideremos em G = (N, E) um separador S que parte N — S com componentes conexas 
Ni,N2, . . . ,Nk. Em cada uma desta componentes podemos considerar um novo separador que 
a reparte, e assim recursivamente. Podemos representar este processo pela árvore de dissecgáo, 
D, onde S é a raiz, uma componente separada por S, ou o separador dentro dela, é filha de S, e 
as componentes náo repartidas sáo as folhas da árvore. 

Uma pós-ordem dos vértices de uma árvore H de raiz r é uma ordem q, dos vértices de H, que 
hsta os vértices de cada uma das árvores de H — r, recursivamente em pós-ordem, e, finalmente, 
r. Seja q uma pós-ordem numa árvore de disseccáo, D, de G. Substituindo em q cada nó, d, de 
D pelos vértices de G em d, temos uma ordem de dissecgáo, q. 



Lema 6.6 (DissecQao) Consideremos a eliminagáo dos vértices de G na ordem de dissecgáo q. 
A eliminagáo de um dado vérice, v E d, só pode preencher lados dentro do seu nó em D, d, ou 
entre (vértices de) d e (vértices em) nós ancestrais de d (em D), ou ainda entre (vértices em) 
ascendente de d (em D). 



Demonstracáo: Trivial, pelo lema do caminho. 

Exemplo 7: 

Vejamos um exemplo de ehminagao usando uma ordem de dissecgáo 

- - - 1 2 {1,2} 

/ / I I \ / I \ 

3 I 4 5 6 7 {3,4} {5,6} {7} 

8-9-10 11 12 13 {8,9,10} {11,12} {13} 
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O último conjunto da partigáo, Sk ou a raiz de G, é em, G, um separador entre cada uma 
das componentes de G — A;. Para minimizar a regiáo de possível preenchimento (em G ou QAQ') 
gostaríamos de ter o separador Sk "pequeno e balanceado" , i.e. tal que 

• "^Sk seja o menor possível. 

• Sub-árvores tenham aproximadamente o mesmo número de vértices de G. 

Recursivamente, gostaríamos de ter como raiz de cada uma das sub-árvores um bom separador, 
i.e., pequeno e balanceado. Veremos a seguir várias heurísticas para obter num grafo qualquer, 
G, um bom separador. 

Heurística de Busca em Largura: 

Uma busca em largura, BEL, a partir uma raiz v & N, particiona os vértices de G em níveis 
Lq, Li, . . . Lfc, definidos por 

Lo = {v} , Li+i = adj{Li) - Li_i . 

A profundidade do nível L^ é i, e & largura do nível Lj é ^Lj. A profundidade e a largura da 
BEL sáo, respectivamente, a máxima profundidade e a máxima largura nos níveis da BEL. 

Lema 6.7 O nível L^ separa, em G, os vértices em níveis mais profundos dos em níveis menos 
profundos que i. 

A heurística de BEL procura um separador balanceado S (^ L^, tomando i ~ k/2, ou entáo 
tomando i \ EÍ"^ #^j < n/2 A E7+i ^Lj < n/2 . 

Para obter um separador pequeno a heurística procura uma raiz, v, que gere uma BEL de 
máxima profundidade, com o intuito de reduzir a largura da BEL. A distáncia (em G) de um 
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vértice v a um vértice w, dist{v, w), é o comprimento, ou número de lados, do caminho mais curto 
entre ambos os vértices. A excentricidade de um vértice v é exc{v) = maxw eNdist{v,w). Um 
vértice de máxima excentricidade se diz periférico, e sua excentricidade é o diámetro de G. 

Uma BEL com raiz v terá profundidade igual a excentricidade de v. Isto motiva querermos 
iniciar a BEL por um vértice periférico. Encontrar um vértice periférico é um problema com- 
putacionalmente difícil. A heurística de Gibbs encontra um vértice quase-periférico como 
segue: 

L Escolha como raiz um vértice de grau mínimo. 

2. Forme os níveis da BEL com raiz v, Li . . .L^. Particione o nível mais profundo em suas 
componentes conexas, L^ = ^{Sj, e tome um vértice de grau mínimo, Vj, em cada compo- 
nente. 

3. Para j = 1 : l 

• Tome Vj como nova raiz e encontre os níveis da BEL, Li . . . L¡¡.i 
Ate que k' > k ou j = l. 

4. Se o passo 3 terminou com k' > k, volte ao passo 2. Caso contrário a atual raiz é um vértice 
quase-periférico. 

Exemplo 8: 

Retomando o exemplo 7, a heuristica de Gibbs encontra 3 como vértice quase-periférico. 
Tomando 3 como raiz geramos a árvore H por BEL: 

10-4 13 

/ / 

H= 3-8-9-1-5-12-6-2-7 

L=12345 6 78 9 

Escolhendo 5*1 = L5 como primeiro separador, e depois S^ = L^ e S^ = L-j como sepa- 
radores dentro de cada uma das componetes separadas por Si, obtemos a odem por disseccáo 
q = [3, 8, 1, 10, 9, 11, 12, 2, 13, 7, 6, 4, 5]. 
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Note que nas linhas (colunas) correspondentes aos vértices do primeiro separador, 5*1 = {1}, 
pode haver ENN's em qualquer posigáo. Note também que o resto da matriz está em forma diag- 
onal blocada (vide definigáo no capítulo 8), onde cada bloco correspnde a uma das componentes 
separadas por 5*1. Esta estrutura se repete em cada bloco, formando a estrutura "espinha de 
peixe" característica de ordens por dissecgáo. Note que esta estrutura é preservada pela fatoragáo 
de Cholesky. 



Exercícios 

1. Implemente a ehminagáo simbóhca num grafo G, com a ordem q, computando F e H, em 
tempo 0{^enn{L)). 

2. Quáo eficiente {k em tempo 0{n^)) poderíamos implementar o algoritmo imphcito na úhima 
parte do teorema de caracterizagáo de grafos cordais? 

3. Considere o seguinte algoritmo para numerar os vértices de um grafo na ordem n,n — 

1,...2,1: 

Ordenamento Reverso por Grau Máximo (ORGM): 

(a) Escolha, como vértice n, um vértice qualquer. 

(b) Escolha, como vértice seguinte um vértice ainda náo numerado adjacente a um número 
máximo de vértices já numerados 

Prove que ORGM define uma ordem perfeita. 

4. Quáo eficientemente podemos implementar o ORGM? 

5. De um exemplo de vértice quase-periférico que náo seja periférico. 
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Capítulo 7 



ESTRUTURA 



Acoplamento de Sub-Sistemas 



Estruturas Blocadas 

Consideraremos neste capítulo matrizes com blocos de elementos nulos dispostos de forma regular. 
Estudaremos duas estruturas: a triangular-blocada superior, e a angular blocada por colunas. 
O bloco IB tem dimensáo m{r) x n{s), e na estrutura triangular blocada m{k) = n{k). 
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Estas estruturas blocadas propiciam grandes facilidades computacionais. Em particular tri- 
angularizá-las corresponde a triangularizar os blocos diagonais, sendo que nenhum elemento náo 
nulo é criado, durante a triangularizacáo, nos blocos nulos. 



7.1 Estrutura Triangular Blocada 

Como já visto no estudo de matrizes de permutagáo, dado G = {N, B), N = {1, . . . n}, B, n x n, 
sua matriz de adjacéncia, e um reordenamento de seus vértices, q = [gi,...g„], g¿ G A^, entáo 
a matriz de adjacéncia do grafo G com os vértices reordenados (i.e., reindexados) por q é B = 

BtS = QBQ'- 

Lema 7.1 Considere o reordenamento coerente dos vértices de G = {N,B) numa ordem coerente 
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Q = [^Qí'^Q ■ ■ -^(l]- Conforme a definigáo de ordem coerente, cada bloco, ^q = ['^qi, ■ ■ -^qnik)], 
contém os vértices de uma CFC (componente fortemente conexa) de G, v^, e {vi,...Vh) estáo 
topologicamente ordenados. Neste caso a matriz de adjacéncia do grafo reordenado, B = QBQ' , 
é triangular hlocada superior, de blocos ^B, n{r) x n{s). 



Uma matriz que náo possa ser, por permutacoes de linhas e de colunas, reduzida, isto é, posta 
na forma triangular-blocada, é dita irredutível. 

Dada uma matriz A, náo singular, procuraremos permutagoes de linhas e colunas, isto é, por 
matrizes de permutagáo R e Q, encontrar A = RAQ que seja a "mais fina" partigáo possível de 
A. Observemos que a hipótese de náo singularidade de A imphca na náo singularidade dos blocos 
diagonais, pois 



det{A) = det{A) = ]J det{lA) . 



k=l 



Sejam R e Q matrizes de permutagao e seja P = QR, 

Á = RAQ = Q-^QRAQ = Q'PAQ 

ou seja, podemos escrever qualquer transformagáo do tipo RAQ como uma permutagáo de hnhas, 
PA, seguida de uma permutagáo simétrica Q'{PA)Q. 

Consideremos agora o grafo associado á matriz PA, G{PA), que tem por matriz de adjacéncia 
o padráo de esparsidade de PA, B{PA), i.e., 

G{PA) = {N, B{PA)) = {N, r), j G Vii) ^ {PA)¡ ^ . 

Do úhimo lema sabemos que a permutagáo simétrica que dá a mais fina partigáo de PA é um 
ordenamento coerente dos vértices de G{PA). Ademais, PA é irredutível se G{PA) é fortemente 
conexo. Resta, portanto, anahsar o papel da permutagáo de hnhas, P, na permutagáo geral 
de hnhas e colunas de A = Q'PAQ, onde a permutagáo simétrica é dada por um ordenamento 
coerente em G{PA). 

Pela náo singularidade, A deve possuir ao menos uma diagonal de elementos náo nulos (náo 
necessariamente a diagonal principal, veja as definigóes de determinante e diagonal). Portanto 
existe uma permutagáo de hnhas, P, que posiciona esta diagonal de elementos náo nulos na 
diagonal principal de PA. Uma tal permutagáo P será dita uma permutaQao própria, ao passo 
que uma permutagáo que coloque algum elemento nulo na diagonal principal será dita imprópria. 
Isto posto, mostraremos que: 
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Teorema 7.1 

1. Qualquer permutagáo própria induz a mesma estrutura de partigáo, isto é, h blocos de di- 
mensóes n(l) . . . n{h), com os mesmos índices em cada bloco. 

2. Qualquer permutagáo imprópria náo pode induzir uma partigáo mais fina em A. 



Demonstragao: 

Pelo lema de Hoffmann G{PA) é fortemente conexo, i.e., PA é irredutível se qualquer conjunto 
de k < n linhas tiver elementos náo nulos em ao menos k + 1 colunas. Como esta caracterizagáo 
independe da ordem das linhas, mostramos que a irredutibihdade da matriz PA independe da 
permutagáo própria considerada. 

Da mesma forma, a irredutibihdade do bloco sudeste, ^A bem como a existéncia da CFC v^ 
entre as "úhimas" (no sentido da ordem natural do grafo reduzido) componentes associadas a 
qualquer outra permutagáo própria, P, está garantida, de modo que a invariáncia da estrutura de 
A segue por indugáo no número de componente fortemente conexas (blocos). 

Mostramos agora que se A tiver algum elemento diagonal nulo, entáo as componentes forte- 
mente conexas de G{A) sáo fusóes de componentes fortemente conexas de G{PA): Distinguamos 
os zeros na diagonal de A e consideremos o grafo G^{PA) obtido de G{PA) ao adicionarmos as 
arestas correspondentes aos zeros distinguidos em PA. As CFCs em G{A) sáo exatamente as 
CFCs em G^{PA), sendo que as arestas adicionais (se quando adicionadas ao grafo reduzido de 
G{PA) formarem ciclos) só podem tornar equivalentes alguns vértices antes em CFCs distintas. 
QED. 

A permutagáo própria, P, pode ser vista como um casamento perfeito entre hnhas e colunas, 
onde casamos a coluna j com a hnha p{j), entre seus pretendentes r^^(j) = {i \ B¡ ^ 0}. 
Podemos portanto encontrá-la através do algoritmo Húngaro, de preferéncia com o emprego de 
alguma heurística eficiente para evitar a freqüente geragáo de árvores de caminhos de aumento. 

Exemplo 1: 

Considere as matrizes de adjacéncia B, sua permutagáo própria PB, e o posterior ordenamento 
coerente Q'PBQ. 
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Q'PAQ = B\ 
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Observe que G{B) é fortemente conexo. Em B destacamos a diagonal a ser posicionada por P 
na diagonal principal de PB, e o zero inicialmente na diagonal principaL Em PB destacamos 
o mesmo zero originalmente na diagonal de B. G{PB) tem 3 CFCs, que vohariam a fundir-se 
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numa só, caso adicionássemos a aresta correspondente ao zero destacado em PB. Neste exemplo 
tivemos p = [3, 1, 2], q = [1, 3, 2]. 

O procedimento P4 (PartÍQao e Pré-Posicionamento de Pivós) explora esparsidade estrutural, 
e posteriormente a esparsidade local a cada bloco, como segue: 



1. Encontre uma permutacao própria, PA, através do algoritmo húngaro complementado com 
uma heurística eficiente. 



2. Encontre, através do algoritmo de Tarjan, um ordenamento coerente Q'PAQ. 



3. Inverta este úhimo ordenamento, QPAQ', de modo a colocar a matriz na forma triangular 
blocada inferior. Em seguida aphque a heurística P3 a cada bloco diagonaL 



Sáo apresentados trés exemplos de aphcacáo do P4. A disposÍQáo original dos ENNs da matriz 
original, A, corresponde aos sinais +, *, ou números de 1 a 9 no corpo da matriz. Como prescrito 
no primeiro passo do P4, primeiramente encontramos uma diagonal náo nula, ou equivalentemente 
uma permutagáo própria PA. O vetor inverso de índices de hnha permutados, p, é dado á esquerda 
da numeragáo original das hnhas. Os asteriscos indicam os elementos desta diagonal. 

Como prescrito no segundo passo do P4, devemos em seguida encontrar um reordenamento 
coerente, q, em G{PA). Para tanto aphcamos o algoritmo de Tarjan: Primeiramente fazemos a 
busca em profundidade canónica em G{PA). Para tanto percorremos as hnhas de PA, gerando a 
primeira fioresta de cada exemplo. As raízes desta busca sáo assinaladas por um acento circunfiexo, 
tanto na fioresta como no vetor p. Os números no corpo da matriz correspondem a ordem de 
visitagáo nesta busca. A inversa da ordem de retorno é dada pelo vetor b. 

O algoritmo de Tarjan exige em seguida a busca em profundidade canonica no grafo inverso 
reordenado por b, G{{BPAB)'); na verdade basta tomarmos as raízes desta busca na ordem 
canonica. Percorrendo as colunas da matriz, construímos a segunda fioresta em cada exemplo, cu- 
jas raízes estáo assinaladas por um acento circunfiexo, tanto na fioresta como no vetor b. A ordem 
de visitagáo nesta segunda busca em profundidade nos dá o reordenamento coerente QPAQ'; ap- 
resentado exphcitamente em cada exemplo para que se reconhega a estrutura triangular superior. 

Finalmente, ao final dos exemplos, apresentamos o terceiro passo do P4: a inversáo deste 
ordenamento coerente, Q'PAQ, (portanto triangular inferior), com posterior aphcagáo do P3 a 
cada bloco diagonaL Neste ponto fica claro que o segundo e o terceiro exemplos diferem apenas 
pela permutagáo original da matriz. Neste ponto indicamos também os zeros preenchidos durante 
a fatoragáo. 
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Exemplo 2: 
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Exemplo 3: 
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Exemplo 4: 
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Exemplos 2 a 4; Passo 3 do P4: 
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7.2 Estrutura Angular Blocada 



Suponhamos dada uma matriz quadrada e náo singular na forma angular blocada, com blocos 
diagonais ^B, . . . ^B, ^B m{k) x n{k), d{k) = m{k) — n{k) > 0, e os correspondentes blocos nas 
colunas residuais ^C . . . ^C , ^C m{k) x n{h + 1). Podemos usar rotacóes de Givens para fatorar 

ky 



cada um dos blocos diagonais, ^B 
bloco de zeros é d{k) x n{k). 



'Q 



onde V é triangular superior n{k) x n{k), e o 



'B 



'C 



'B ^C 





'V 




'w 









'z 


5 




hy 


^w 









^z 



Para completar a fatoragáo QR da matriz original, respeitando a estrutura de blocos, permu- 
tamos os blocos ^Z para as últimas linhas da matriz, formando o bloco quadrado Z de dimensáo 
Y,id{k) = n{h + 1). Finalmente completamos a fatoracáo QR do bloco sudeste, Z = QS. 



'V 



^W 



hy hyy 

'z 



'V 



^w 

hy hyy 

s 



Como permutacóes sáo apenas um tipo especial de transformacóes ortogonais, obtivemos o 
fator triangular da fatoragáo QR respeitando a estrutura diagonal blocada da matriz original, 
B = QU . A inversa da matriz original seria dada por B^^ = U^^Q', mas usando que Q' = U^^B' 
temos B^^ = U^^U^^B' . Isto é, obtivemos uma fatoracáo de B onde todos os fatores herdam (e 
sáo computados de acordo com) a estrutura angular blocada da matriz original. 

Observando que B'B = U'Q'QU = U'U , vemos que uma maneira alternativa de computar o 

fator triangular da fatoragáo ortogonal de B, é computar o fator de Cholesky da matriz simetrizada 

B'B: 

^B'^B ^B'^C 



ir-'/i 



C'^B 



hQihQ 

h(J,hQ 



h^ihfj 



Ao ehminarmos os h blocos das hnhas residuais de B' B, formamos o bloco sudeste Z = Z+^^ Z, 
a ser fatorado na úhima etapa do processo, Z = S'S. Ao final, obtemos exatamente o fator 
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triangular da fatoracáo QR da matriz original. 



7.3 Particáo de Hipergrafos 

Um hipergrafo é um par ordenado G = {V, C) onde o primeiro elemento é um conjunto finito, 
o conjunto de vértices, e o segundo elemento é uma matriz booleana, a matriz de incidéncia. 
Se \V\ = m, cada coluna de C, m x n, nos dá os vértices sobre os quais incide o correspondente 
(hiper)lado. No caso particular de todos as colunas terem exatamente 2 ENN's temos na matriz 
de incidéncia mais uma representagáo de um grafo simétrico. Em hipergrafos todavia um lado 
genérico pode incidir sobre mais de 2 vértices. 

O problema de permutar PAQ para forma angular blocada pode ser visto como um problema 
de partigáo no hipergrafo G = {M,B{A)), M = {l,...m}, B{A) a matriz booleana associada 
á matriz A, m x n. No problema de partigáo damos a cada hnha i G M de B{A) uma cor 
p{i) & H = {1, . . . h}. A cor de cada lado é definida como o conjunto de cores dos vértices sobre os 
quais este incide, q{j) = {p{i) \ Aj ^ 0}. Lados muhicoloridos correspondem a colunas residuais 
na forma angular blocada, e os lados de cor q{j) = k correspondem ás colunas no bloco angular 
formado pelos ENN's Aj \ p{i) = k A q{j) = {k}, k e H = {1, . . .h}. 

Para definir o problema de partigáo faha-nos uma fungáo objetivo a ser minimizada. Em vista 
das aphcagóes já apresentadas, e outras a serem apresentadas no próximo capitulo, queremos ter: 

• Aproximadamente o mesmo número de linhas em cada bloco. 

• Poucas colunas residuais. 

Com estes propósitos é natural considerar a seguinte fungáo de custo de uma dada partigáo de 
hnhas em cores p : M ^ H: 

h 

f{p) = c{p) + a^{m/h - s{k)f , 

k=l 

onde c{p) = \{j G A^ | \q{j)\ > 2}| ,e s{k) = \{i e M \ p{i) = k}\ . 

Mesmo casos especiais deste problema sáo NP-difíceis. Por exemplo: Seja A a matriz de 
incidéncia de um grafo, m um múhiplo exato de /i = 2, e faca a suficientemente grande para 
garantir que todos os blocos tenham exatamente m/h linhas. Este é o problema exato de 2- 
partigáo em grafos, e a versáo de reconhecimento deste problema é NP-Completa; veja problema 
ND14 em [Garey79]. Um algoritmo de anulamento simulado com perturbagóes métricas para 
resolver este problema é apresentado em [Stern92]. 
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7.4 Paralelismo 

Um dos fatores mais importantes no desenvolvimento de algoritmos é a possibilidade de realizar 
várias etapas de um procedimento em paralelo. No restante deste capítulo adaptaremos algumas 
das fatoragóes anteriormente estudadas para as estruturas blocadas, visando paralelizar etapas 
independentes. Vejamos a seguir alguns conceitos básicos de computacáo paralela. 

Existem vários modelos teóricos de computador paralelo, e inúmeras instáncias e implementagóes 
destes modelos em máquinas reais. O modelo mais simples é o de memória compartilhada. 
Neste modelo vários processadores, tém acesso a uma memória comum. Neste modelo, a descrigáo 
de uma algoritmo paralelo envolve basicamente dois fatores: 

• Como distribuir o trabalho entre os processadores. 

• Como sincronizar as diversas etapas do algoritmo. 

Este modelo é conceitualmente simples e elegante; todavia limitagoes da nossa tecnologia inviabi- 
lizam a construgáo de máquinas de memória compartilhada com mais de uns poucos (da ordem 
de dez) processadores. 

O modelo de rede é mais genérico. Nele o computador é visto como um grafo: cada vértice, 
nó, ou processador representa: um processador propriamente dito, uma memória local, i.e., 
acessível somente a este processador, e portas de comunicagáo. Cada aresta representa uma via 
de comunicagáo inter-nós. Note que no modelo de memória compartilhada, a comunicagáo entre 
os processadores podia ser feita de maneira trivial através da memória; todavia no modelo de rede 
é preciso saber os detalhes da arquitetura da máquina para especificar um terceiro aspecto do 
algoritmo: 

• A comunicagáo entre os processadores. 

Estes detalhes incluem a disposigáo das vias de comunicagáo, ou topologia, a velocidade de 
comunicagáo em relagáo a velocidade de processamento, a possibilidade ou náo de haver comu- 
nicagoes simuháneas em vias distintas, etc. Algumas destas topologias comumente empregadas, 
sáo: Estrela, Barra, Anel, Grade, Toro, Hipercubo e Borboleta. 

Como exemplo de algoritmo paralelo, calculemos a média de n números numa rede com p 
processadores. Suponhamos que inicialmente tenhamos n/p destes números em cada uma das 
memórias locais. Por simplicidade suponhamos que n é um múhiplo de p, e que n » p. Na 
primeira fase do algoritmo cada processador, k, calcula a média dos n/p números em sua memória 
local, m{p). Se os processadores sáo todos iguais (rede homogénea), cada processador completa 
sua tarefa em 1/p do tempo necessário para calcular a média geral num computador com apenas 
um processador deste mesmo tipo. Na segunda fase do algoritmo reunimos as médias parciais para 
calcular a média geral, m(0) = (l/p) YTk=i^{k)- Examinemos como calcular esta média geral em 
duas redes com topologia de anel, onde cada qual: 
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1. Náo permite comunicagóes em paralelo. 

2. Permite comunicagóes em paralelo via segmentos de arco náo superpostos. 

Novamente por simplicidade, suporemos que p = 2^^. 
Na primeira rede, para k = 1 : p, 

1. Calcule no, nó k, s{k) = s{k — 1) + m{k). 

2. Transmita s{k) ao processador k + 1. 

Neste procedimento temos as somas parciais das médias s{k) = J2'í=im{k), a condigáo de inicial- 
izagáo é s(0) = 0, e ao término do algoritmo podemos computar, no nó p, a média m(0) = s{p)/p. 

Na segunda rede, para i = 1 : q, 

• j = 2*; em paralelo, para k = j : p, 

1. Calcule, no nó k, r{i, k) = r{i — l,k) + r{i — 1, /c — j/2). 

2. Transmita r{i, k) do nó k para o nó /c + j. 

Neste procedimento temos as somas parciais das médias r{i,k) = Y!¡=k-jj^im{l), a condigáo de 
inicializagáo é r(0, /c) = m{k), e ao término do algoritmo podemos computar, no nó p, a média 
m(0) = r{q,p)/p. 

Em geral, a transmissáo de dados entre nós é muito mais lenta que a manipulagao destes dados 
localmente e, ao medir a complexidade de um algoritmo, contamos separadamente os trabalhos 
de processamento e comunicagáo. 



7.5 Fatoragoes Blocadas 

Das segóes anteriores vemos que quase todo o trabalho na fatoragáo QR de uma matriz angular 
blocada, A = QU, ou da fatoragáo de Cholesky A'A = U'U , consiste na aphcagáo repetitiva de 
algumas operagóes simples sobre os blocos. Para tirar vantagem desta modularidade em algoritmos 
para fatoragáo e atuahzagáo de matrizes com estrutura angular blocada, definimos a seguir algumas 
destas operagóes, e damos sua complexidade em número de operagóes de ponto fiutuante. 

1. Compute a fatoragáo de Cholesky parcial, ehminando as primeiras n colunas da matriz 
blocada 

' F G 

G* 
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para obter 



V W 
Z 



onde F = F^énxn, eGénxl. Isto requer (l/6)n^ + (l/2)n^/ + (l/2)n/^ + O^n"^ + P) 
FLOPs. 

2. Compute a transformagáo inversa parcial, i.e. m, em 





t 


' u' ' 

_u\ 


= 


_y\ 



onde K é n X n triangular superior, W^énx/, Oe/ sao as matrizes zero e identidade, e u 
e y sáo vetores coluna. Isto requer (l/2)n^ + n/ + 0(n + /) FLOPs. 

3. Reduza a triangular superior \xma, matriz de Hessenberg, i.e., aplique a seqüéncia de rotagóes 
de Givens G(l, 2, 9), G{2, 3,6) .. . G{n — 1, n, 9) á matriz blocada 

V W 

onde V é nxn Hessenberg superior, eW é nxl, para reduzir V a triangular superior. Isto 
requer 2n^ + Anl + 0{n^ + f) FLOPs. 

4. Reduza a triangular superior uma matriz blocada coluna - triángulo superior, i.e., aplique 
a seqüéncia de rotagóes de Givens G{n — 1, n, 9), G{n — 2,n — 1,9), . . . G{1, 2, 9) á matriz 
blocada 

' u V' 

onde u é um vetor coluna nxl,el^énxn triangular superior, de modo a reduzir u á 
um único ENN na primeira linha, assim transformando V de triangular para Hessenberg 
superior. Isto requer 2n^ + 0{n) FLOPs. 

Em ambas as fatoragóes, A = QU e A'A = U'U , muitas das operagoes nos blocos podem ser 
feitas independentemente. Portanto a estrutura angular blocada náo só nos dá a possibilidade de 
preservar esparsidade, mas também a oportunidade de fazer várias operagóes em paralelo. 

Descreveremos uma forma de paralelizar a fatoragáo de Cholesky numa rede de /i+ 1 nós. Para 
k = 1 . . .h alocamos blocos das matrizes A e U a. nós específicos, como segue: 

• Os blocos D^ E^, V^ e W^ sáo alocados ao nó k. 

• Os blocos sudeste, Z e S, sáo alocados ao nó (ou h + 1). 



Expressaremos a complexidade do algoritmo em termos da soma e do máximo das dimensoes 
dos blocos. 



dhsum = y^^m{k) 
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dbmax = max{m(l), . . . , m{h),n{h + 1)}. 

Na análise de complexidade contabilizaremos o tempo de processamento, medido em FLOPs, 
pTime, bem como a comunicagáo inter-nós, INC. Quando h operagóes sobre blocos, bop^ . . . bop^, 
podem ser efetuadas em paralelo (em nós distintos), contamos seu tempo de processamento por 
f\^flops{bop^) = flops{bop^) A . . . A flops{bop^) , onde A é o operador máximo, e flops{bop'') 
é o numero de operagóes de ponto flutuante necessário para a operagáo no bloco k, bop{k). Nas 
equagóes que seguem, A tem precedéncia menor que qualquer operador multiplicativo ou aditivo. 
As expressóes "No nó k=l:h compute" ou "Do nó k=l:h envie" significam, "Em (de) todos os nós 
1 < k < h, em paralelo, compute (envie)". Nas expressóes de complexidade ignoraremos termos 
de ordem inferior. 

Damos agora uma descrigáo algorítmica da fatoragáo de Cholesky blocada bch{): 

1. No nó k=l:h compute os blocos {B^fB'', {B^fC^, e {C'^fC^. 
pTime = m{k)n{k)'^ + m{k)n{k)n{h + 1) + m{k)n{h + 1)^ < Sdbmax^, 
INC = 0. 

2. Envie (C'^')*C'^ do nó k para o nó 0, onde acumulamos Z^ = Y^\{C^YC^. pTime = h n{h + 
1)2 < h dbmax"^ , INC = h n{h + 1)^ < h dbmax"^ 

3. No nó k compute a fatoragao de Cholesky parcial, ehminando as primeiras n{k) colunas, da 
matriz blocada 

{B'^yB^ {B^yc^ 

{C^'B^ 

obtendo 

yk y^k 

z^ 

pTime = {l/Q)n{k)^ + {l/2)n{k)'^n{h + 1) + {l/2)n{k)n{h + 1)^ < {7/6)dbmax^, 
INC = 0. 

4. Envie Z'^ do nó k para o nó 0, onde acumulamos Z = J2o Z'^. 
pTime = h n{h + 1)^ < /i dbmax"^ , INC = h n{h + 1)"^ < h dbmax"^ . 

5. No nó fatore o bloco sudeste S = chol{Z), onde chol{) indica a fatoragáo de Cholesky 
padráo. 

pTime = {l/6)n{h + 1)^ < {l/6)dbmax^ , INC = 0. 

Teorema 7.2 A fatoragáo de Cholesky blocada, bch{), requer náo mais de (4 + l/3)dbmax^ + 
h dbmax"^ tempo de processamento, e h dbmax"^ tempo de comunicagáo inter-nós. 

Nos passos 2 e 4, se a rede permite comunicagóes em paralelo, a reuniáo da matriz acumulada 
pode ser feita em log{h) passos, e podemos substituir h por log{h) no úhimo teorema. 



Capítulo 8 



ESCALAMENTO 



Representagáo em Ponto Flutuante 



8.1 O Sistema de Ponto Flutuante 

A representaQáo de um número real, ( & R, em um computador tem, usualmente, precisáo finita. 
A inexatidáo desta representagáo introduz erros no resultado final do processamento e o objetivo 
desta segáo é obter limites máximos para estes erros quando da aplicacáo do método de Gauss. A 
representagáo normalmente utilizada para números reais é o sistema de representagáo em ponto 
flutuante normalizado de t dígitos e base b, SPF, isto é, 

//(C) = ±0.diá2...áí*&^",OU 

±0.did2 . . .dt E ±n 



onde 



dk,n, b E N 
0<dk<b, diy^O 
< n < emax, b j^ . 



Dado um real z/ com representagao exata num dado SPF, a fragao normalizada será denominada 
mantissa. A mantissa e o expoente de z/ seráo denotados, respectivamente, 

mant{v) = ±0.di . . .dt, 
expo{v) = ±n 

Há duas maneiras normalmente empregados para obter fl{C) a partir de (: truncamento 
e arredondamento. Para trunc{() simplesmente truncamos a representagáo em base b de ( 
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obtendo uma mantissa de t digitos. No arredondamento tomamos a mantissa de t digitos que 
melhor aproxima (. Assim, para C = n = 3.141592653..., e o SPF com 6 = 10 e í = 5, 
trunc{C) = +0.31415 E + le round{C¡ = +0.31416 E + 1. 

Note que C = ^^o pode ser representado devido a condigáo di j^ 0. Portanto, alguma 
representagáo inambígua deve ser convencionada, por exemplo, //(0) = +0.00 . . . £" + 0. Se n > 
emax náo há representacáo possível no SPF, e dizemos que houve "overflow" ou transbordamento. 

8.2 Erros no Produto Escalar 

Definimos a unidade de erro do SPF, u, como 6^~* no caso de usarmos truncamento, e b^~^/2 
no caso de usarmos arredondamento. 

Lema 8.1 Dado ( E R e fl{() sua representagáo, num dado SPF de unidade de erro u, 35 G 

[-uM |//(c) = c(i + 5). 



Demonstracao: 

Pela definigáo de SPF, se expo{fl{()) 

mo-ci 

ICI 



e, ve-se que 



\6\< 



uh' 



e-l 



Ue~l 



U 



Se z/ e /i sáo números em ponto flutuante, isto é números reais com representagáo exata num 
dado SPF, é perfeitamente possível que uma operagáo aritmética elementar entre eles resulte num 
número sem representagáo exata. Do lema anterior, porém, sabemos que, qualquer que seja a 
operagáo aritmética, -k G {+, — , *, /}, fl{u */i) = (z/t<c/í)(1 + 5), para algum 5 | |5| < u. QED. 

Exemplo 1: 

Consideremos um computador que armazena um número real em 4 bytes, sendo 3 bytes para 
os dígitos da mantissa e 1 byte para o sinal do número, o sinal do expoente e os 6 bits restantes 
para o módulo do expoente como um inteiro em base 2. Supondo que todos os cálculos sáo feitos 
com arredondamento, calculamos, a unidade de erro do SPF, -u, e o maior real representável, 
rmax = b^'^"'^, no caso de usarmos base b = 256, 6 = 16 ou 6 = 2. 



b 


t 


u = b'~'/2 


rmax = b^'^ 


2 


24 


2-24 < iQ~7 


264 > ^^19 


16 = 2^ 


6 


2-21 < 10-6 


2256 > 10^5 


256 = 2^ 


3 


2-17 < 10-5 


2512 > iQlSO 



Um recurso freqüentemente disponível, concomitantemente ao uso de um SPF de base b e t 
dígitos, é o sistema de representagáo em ponto flutuante normahzado de precisáo dupla, SPFD, 

com mantissa de 2í dígitos, que denotamos fld{(). 
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Este recurso é extremamente útil, se utilizado parcimoniosamente. Podemos trabalhar com a 
maior parte dos dados no SPF, de precisáo simples, e utilizar a precisáo dupla apenas nas passagens 
mais críticas do procedimento. Em analogia ao SPF, definimos a unidade de erro do SPFD, ud 
como 6^~^* no caso de usarmos truncamento, e 6^~^*/2 no caso de usarmos arredondamento. 

É interessante notar que se z/ e /i sáo números reais com representagáo exata no SFP, de precisáo 
simples, a representagáo do produto destes números no SPFD é exata, isto é fld^u * fi) = u * fj,, 
pois o produto de dois inteiros de t dígitos tem, no máximo, 2í dígitos. 

Estudamos agora o efeito cumulativo dos erros de representagáo, em todas as passagens in- 
termediárias num produto interno. Sejam x, 1 x n e y, n x 1, vetores cujas componentes tém 
representagáo exata, num dado SPF. Definimos 

flixy) = fliflixnyn + fliflixn-iy""') + ... + fliflix^y') + flixiy')) ...)). 
Lema 8.2 Dados u & [0, 1[, n E ISí tq nu < a < 1 , e Si \ |5¿| <u, i = 1 . . .n, entáo 

n 

l-nu< [|(1 + Si) <l + il + a)nu . 
1 

Demonstragáo: 
Como 

n 

(l-nM)"<n(l + í^)< (l + nw)", 

1 

basta provar que 

• (1 — m)" > 1 — nu. 

• (1 + m)" < 1 + nu + anu. 

Considerando a fungáo /(-u) = (1 — «)"" temos que 36 g]0, 1[ | 

fiu) = fi0) + uf'iu)+u'f"ieu)/2 

= 1 - nw + n(n - 1)(1 -^m)""V/2 . 

Da náo negatividade do último termo segue a primeira inequagáo. 
Considerando que V/3 G [0,a], 

l + (3<e^<l + (3 + a(3, 

temos que 

(1 + -u)" < emu) < 1 + nw + anu . 

QED. 
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Corolário: 

Nas condigoes do lema 2, 39 E [—1, 1] tq 

n 

1[{1 + 6i) = 1 + 9(1 + a)nu . 
1 

Lema 8.3 Se x, 1 x n e y, n x 1, sáo vetores cujas componentes tem representagáo exata num 
dado SPF, de unidade de erro u, com nu < 1, entáo 

n 

fl{xy) =Y,Xiy\l + 0i{l + a){n-i + 2)u) . 

i=l 

Demonstracáo: 

36i, \6i\ < ud, e 6i, \0i\ < 1 \ 

fl{xy) = fl{fl{xr,y-) + fl{fl{x^.w''-') + ... 

+{xsy'{l + 6,) + {x^y^l + 6^) + x,y\l + 6,)){1 + 6,)){^ + 6,).. .)) 
= fl{fl{xnyn + fl{fl{xn-iy^-') + ■■■+ x,y\l + 6,){1 + 6,) 

+W(1 + 52)(1 + 6s){l + 6,) + x,y\l + 6,){1 + 6s){l + 6^) ■ ■ ■)) 

= X„í/"(l + (52„„2)(l + 52n-l) + 

+X„_1Í/"~^(1 + 52n-4)(l + 52n-3)(l + ^^n-l) + • • • 
+X2y\l + (52) (1 + (53) ... (1 + (52„-3)(l + 62n~l) 
+Xiy\l + 5i)(l + 53) ■ ■ ■ (1 + 52„-3)(l + 62n-l) 

= x„y"(l + 9r,{l + a)2u) + a;„„iy"-i(l + í^„_i(l + a)3u) + ... 
+X2y^{l + ^2(1 + a)nu) + xiy\l + ^1(1 + a){n + l)u) . 

Lema 8.4 Se x, 1 x n e y, n x 1, sáo vetores cujas componentes tém representagáo exata num 
SPF, em precisáo simples, de unidade de erro u, sendo n * ud < '-/ < 1, entáo 

fld{xy) = fld{fld{xr.yn + fld{fld{x^^,y^~^) + ... + fld{fld{x2y^) + fld{x,y')) . . .)) 

n 

= ^a;¿í/^(l + ^¿(l+7)(n-¿ + l)Má) 
1 

Demonstracáo: 

3e¿, |e¿| < ud, e 9^, \d.¡\ < 1 \ 

fld{xy) = fld^Xny^" + //á(a;„_ií/"~i + . . . 

+ (x3í/^ + {X2y^ + xiy^){l + ei))(l + 62) . . .)) 
= a;„?/"(l + e„_i) + Xn-iy'''^{l + e„_2)(l + e„_i) + . . . 

W(l + ei) . . . (1 + e„_i) + x^y\l + e^) . . . (1 + e„_i) 
= x„?/"(l + en{l + l)ud) + Xn-iy'''\l + 9n-i{l + i)2ud) + 

X2y^{l + ^2(1 + 7)(^ - ^)ud) + xiy\l + ^1(1 + -f)nud) 
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É freqüente o cálculo em dupla precisao, de produtos de vetores armazenados em precisao sim- 
ples, e o posterior armazenamento do resultado em precisáo simples, isto é o cálculo fl{fld{xy)). 

Do último lema vemos que 

fl{fld{xy)) = {1 + 6)J2 ^.1/Xl + ^.(1 + l){n-t + l)ud) . 

1 

Observagáo 8.1 Se xy ^ J2xiy^^i{^ + l){n — i + l)ud, o que ocorre se n-u ^ 1 e náo houver 
"cancelamentos críticos" na somatória, o resultado final do produto é afetado de um erro da ordem 
do erro introduzido por uma única operagáo aritmética, em precisáo simples. Assumiremos esta 
hipótese no restante do capítulo. 

Exemplo 2: 

Calculemos, no SPF decimal de 2 dígitos com arredondamento, fl{xy), fld{xy) e fl{fld{xy)), 
onde X = [7.5, 6.9, 1.3] ey = [0.38, -0.41, 0.011]'. 

fld{xy) = fld{2.85 - 2.829 + 0.0143) = //á(0.021 + 0.0143) = 0.0353 . 

fl{fld{xy)) = //(0.0353) = 0.035 . 

fl{xy) = fl{2.9 - 2.8 - 0.014) = //(0.1 + 0.014) = 0.11 . 



8.3 Escalamento de Matrizes 

Da anáhse de erros no produto interno, e considerando que a maioria das operagóes nos algoritmos 
de triangularizagáo podem ser agrupadas na forma de produtos internos, fica evidente que é 
conveniente termos todos os elementos da matriz da mesma ordem de grandeza, i.e., termos a 
matriz bem equilibrada. Evitamos assim ter muhiphcadores muito pequenos e soma de termos 
de ordem de grandeza muito diferentes. Se tal náo ocorre para a matriz de um dado sistema, 
Ax = b, podemos procurar x através da solugáo de um outro sistema melhor equihbrado. 

Se E = diag{ei, e^, . . . e„) e D = diag{di, d^, . . . dn), onde ej, di j¿ 0, o sistema EADy = Eb 
é obtido muhiphcando-se a i-ésima equagáo do sistema por ej, e efetuando-se a substituigáo de 
variáveis x = Dy, o que equivale a muhiphcar a j-ésima coluna de A por dj. Uma transforamagáo 
A = EAD, A m X n, E e D diagonais com ej, di ^ 0, é um escalamento da matriz A. 

Exemplo 3: 



EAD 



1 " 




' 1 1 1 ' 




2 




1 1 1 




3 




111 





11 













12 





= 








13 





11 


12 


13 " 


22 


24 


26 


31 


36 


39 
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Estudaremos o problema de escolher um escalamento que "melhor equihbre" uma dada matriz. 
Em primeiro lugar, vale notar que, estando num SPF de base b, a escolha de E e D da forma 
E¡ = lf\ Dj = b'^^ , onde e e d sáo vetores de elementos inteiros, é muito conveniente, pois 
A = EAD e A teráo elementos de mesma mantissa, sendo o efeito do escalamento apenas o de 
aherar os expoentes dos elementos da matriz, 

mant{Al) = mant{Al) , expo{Al) = expo{Al) + e^ + d^ . 

Observagáo 8.2 Para implementar eficientemente a funcáo expoQ e a operacáo de soma de 
um inteiro ao expoente de um número real, devemos conhecer detalhadamente o SPF usado e 
manipular diretamente os campos de bits envolvidos. 

Uma maneira de medir o grau de desequihbrio de uma matriz, A é atravez da média e da 
variáncia dos expoentes de seus elementos: 

mex{A) = ^ expo{Al) / enn{A) , 

vex{A) = ^ (expo{Al) — mex)'^ / enn{A) . 

Nas somatórias que definem mex e vex excluimos os elementos nulos da matriz, os elementos nulos 
da matriz podem ser ehminados das operagóes de produto interno. 

Para náo sobrecarregar a notagáo, doravante escrevemos 

m n m,n 

E'= E > E'= E > E'= E • 

« ¿=1 \A¡yíO J j=l \A¡jíO *'i i,j=l \A\+0 

O primeiro método de escalamento que estudaremos é justamente o método da redugáo de 
variáncia em que procuramos minimizar a variáncia dos expoentes de EAD. 

Tomemos as matrizes de escalamento esquerda e direita como, respectivamente, E = diag{round{e*)) 
e D = diag{round{d*)) , sendo os vetores e* e d* argumentos que minimizáo a variáncia dos ex- 
poentes da matriz escalada. 



vex 



{EAD) = J2 '{expo{AÍ) + ei + dj - mex(A))^ . 



«j 



Um ponto de mínimo deve obedecer ao sistema 
dvex{EAD) 



de* 
dvex{EAD) 



= 2 ^ '{expo{A¡) + e* + d* - mex{A)) 



dd* 



= 2 ^ '{expo{A¡) + e* + d* - mex{A)) 
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ou, fazendo a substituicao ef = e* — mex{A)/2, df = d* — mex{A)/2 

E'et + dt = -j:'expoiA¡) 



J2'4 + d¡ = -J2'expoiA¡ 



Se a matriz A nao tiver elementos nulos, entao a solugao e*, d* do sistema é imediata: 

e* = (1/n) y^(mex(A) — expo{A-¡)) 

j 
d* = (1/m) y^{mex{A) — expo{Al)) 

i 

Esta equagáo pode nos dar uma boa aproximagáo da solugáo em matrizes densas, i.é, com 
poucos elementos nulos, como se considerássemos o "expoente dos elementos nulos" como sendo 
expoente médio de A. Este é o método aproximado da redugáo de variáncia. Procuremos 
agora métodos heurísticos para determinagáo de escalamentos, que sejam menos trabalhosos que 
o método da redugáo de variáncia. 

O método da média geométrica consiste em tomar 

dj = — m¿((y^ ' expo{A¡)) / enn{A^)) 

i 

a = -int{{J2'{expo{AÍ) + dj))/enn{A')) 
j 

Assim os fatores de quihbramento sáo uma aproximagáo do inverso da média geométrica dos 
elementos náo nulos das colunas, ou hnhas, i.e., 

dj = -tnt{{J2"íog{\A¡\))/enn{A^)) 

i 

e, = -tnt{{J2'{\og{\A¡\) + d,))/enn{A')) 
j 

Uma variante do método da média geométrica é o método da média max-min, no qual 
tomamos 

dj = —int{{ma:x'expo{A¡) + min'expo{A¡))/2) 

i i 

e¿ = —int{{max' {expo{A¡) + dj) + min' {expo{A¡) + dj))/2) 
j j 

Finalmente, o método da norma infinito consiste em tomar 

ej = — max' expo{A¡) 
j 

dj = — max' {expo{A¡) + ei) 
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A escolha do particular método a ser empregado depende bastante do tipo de matriz a ser 
equilibrada e da exigéncia que temos sobre vex{A). Em pacotes comerciais de otimizagáo é comum 
a aplicagáo do método max-min um número pré-determinado de vezes, ou até que variáncia dos 
expoentes se reduza a um valor limite aceitável. Este limite deve ser tomado em fungao das 
condigóes do problema, como por exemplo o número de condigáo da matriz, ser definido no capítulo 
9, e da unidade de erro do SPF. 

Exemplo 4: 

Equilibremos a matriz A, dada abaixo, num SPF de base 10, 

1. pelo método aproximado de redugáo de variáncia, 

2. pelo método da média geométrica, 

3. pelo método max-min, 

4. pelo método da norma oo. 

Para A e para cada um dos escalamentos, calculemos mex, vex, e o diámetro da matriz, definido 
como a diferenga entre o maior e o menor expoente dos elementos de A. 



A 



1 


0.7E- 


-4 


0.5E - 1 0.9E0 


-0.2E2 








0.3E-1 


0.3EA 


1 


-0.8E-3 





0.3E6 





0.3E-1 


-0.8E0 0.1E3 


0.7 E9 


:po 


iA) = 


" 

X 



X 


-4-102' 
X -1 X A 
—3 X X 6 
-10 3 9 




+1 
160 
13 



Pelo método aproximado de redugao de variáncia, temos e, d, expo{EAD) e [mex^vex^diam]. 



respectivamente 



mí(8/5) = 2 

mí(2/5) = 

mí(2/5) = 

mí(-6/5) = -1 



int{A/4:) = 1 
mt(12/4) = 3 
mí(6/4) = 2 
mí(l/4) = 
mt(-17/4) = -4 



3 1 


3 


2 





1.57 


X X 




1 

X 


X 
X 



2 


23.4 

4 


X 3 


1 


2 


4 



Analogamente, pelo método da média geométrica, temos 



-mt(4/5) = 1 
-mí(-l/2) = 1 

-mí(l/3) = 
-mt(8/4) = -2 



-mí(0/2) = 
-mí(-8/3) = 3 
-mí(-2/4) = 1 
-mí(3/2) = -2 
-mt(21/4) = -5 



10 1 1-2 

X X 1 X 

X a; 1 

X -1 -1 2 



0.214 

14.4 

4 



8.3. ESCALAMENTO DE MATRIZES 



93 



Analogamente, pelo método max-min, temos 



-mí(-4/2) 
-mí(-2/2) 
-mí(0/2) = 
-mí(4/2) = 



2 

1 



-mí(0/2) = 
-mí(-5/2) = 3 
-mí(-l/2) = 1 
-mí(3/2) = -2 
-mí(ll/2) = -6 



2 1 


2 




-2 


0.143 


X X 




1 

X 


X 
X 


-1 



14.3 
4 


X 


-1 




1 



Analogamente, pelo método da norma-oo, temos 




1 







1 



-3 
-9 






-3 


-1 


-3 


—7 


-1.64 


X 


X 





X 


-4 


59.2 





-2 


X 


X 


-3 


7 


X 















Escalamentos visando equilibrar as matrizes do problema sáo uma etapa importante na solugáo 
de sistemas lineares de grande porte. Note que operagóes de escalamento náo afetam os elementos 
nulos de uma matriz, e portanto náo alteram sua estrutura e esparsidade. O desempenho das 
heurísticas estudadas variam conforme a área de aphcagáo; vale pois testar experimentalmente as 
heurístcas escalamento e suas variagóes. 
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Capítulo 9 



ESTABILIDADE 



Erros e Perturbagoes da Solugáo 



9.1 Normas e Condigáo 

Uma norma, num dado espago vetorial E, é uma fungáo 

1 1 . 1 1 : E ^R\\/x,y e E,a eR : 

1. ||x|| > 0, e ||a;|| = <í^ x = 0. 

2. ||aa;|| = \a\ \\x\\. 

3. ||a; + í/|| < ||a;|| + ||?/||, a desigualdade triangular. 

Sáo de grande interesse em R" as j>-normas 



n 

\x\\„={J2\x,\P^'/P 



\p \ / j 1"^* I 
1 



e de particular interesse a norma 1, a norma 2 (ou norma quadrática, ou Euclidiana) e a norma 
j9 = +00. No caso da norma infinito, devemos tomar o limite da definigáo de j9-norma para 

•p — 7- +00, ou seja, 

I I I I ""11 
X 71 — liia-A. X ''i . 

II IIP ¿=1 I *i 

Dado um espago vetorial normado (E, \\ ||) definimos a norma induzida sobre as trans- 
formagóes lineares limitadas, T : i? — )■ i? tq 3a G R | Va; G -E, ||T(a;)|| < a\\x\\ como sendo 

||T|| = maxí ||T(a;)|| / ||a;|| ) 

II II /.^VII \ /11/ II 11/ 
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ou equivalentemente, por linearidade 

||T|| = max ||T(a;)|| . 
x\ |N|=i" 

Em (R", II II) falamos da norma induzida sobre as matrizes n x n como sendo a norma da 
transformagáo associada, isto é \\A\\ = \\T\\, onde T{x) = Ax, 

Lema 9.1 A norma induzida sobre as matrizes em (R", || ||), goza das propriedades, para 
VA, B n X n, o; G R, 

1. \\A\\ >0 e \\A\\ = 0^A = 

2. \\A + B\\ < \\A\\ + \\B\\ 

3. \\AB\\ < \\A\\ \\B\\ 

Lema 9.2 Para a norma 1 e para norma oo temos as seguintes expressóes explícitas da norma 
induzida sobre as transformagóes, ou matrizes. 



\A\\i = max¿|A|| 



¿=1 

n 



Plloo = maxJ2\A¡\ 

Demonstracáo: 

Para verificar a validade das expressóes observe que 

n n n n 



\\Ax\\, = ElE^i^. l<EEl^il 



\x 



J\ 



¿=1 j=l i=l j=l 



n n 

n 



< Eki|maxEl^ÍI = Plli Iklli 



3=1 ^ ^ i=l 



||^a;||oo = max | y^ A¿Xj |< maxy^ \Al\ 
j=i j=i 



\Xj\ 



n 
n I I n \-^ 1^71 II II 1 1 ^ 1 1 

< max IXjI max > |/1¿ I = ||a;||oo ll^lloo 

7 = 1 i=l ^-^ 

J = l 



e que, se /c é o índice que realiza o máximo na definigao da norma, entao as igualdades sao 
realizadas pelos vetores x = I^ , para a norma 1, e x \ Xj = sig{Al), para a norma oo. 

Definimos o número de condigáo de uma matriz como 



cond{A) = \\A\\ \\A 



-ii 
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Exemplo 1: 

Calcule a norma da matriz de binomial de dimensáo 3, e de sua inversa, nas normas 1 e 
oo. Para cada uma das normas calculadas exiba um vetor x que, multiplicado pela matriz, seja 
"esticado" de um fator igual á própria norma, i.e. x \ \\Ax\\ = ||A||||x||. 



'B 



1 

2 

3 




1 
1 
2 
1 
3 
1 



'B 



' 1 


1 


" 


1 


2 


1 


1 


3 


3 



'B-' 





2 
2 
3 
2 



3 -3 1 
-2 3 -1 

1 -2 1 



Para ^B a soma da norma dos elementos \Al\ por linha e por coluna sáo, respectivamente. 



7 6 4 J e [ 6 8 3 
para a norma 1, ||^-B||, 
norma oo, temos 7, 7, e 49. Finalmente 
rados. 



. Para '^5 ^, analogamente, temos 



7 6 4 



7 6 4 



. Assim, 



B ^ll, e condi^B) sao respectivamente 6, 8, e 48. Analogamente, para 

1 ' r 1 ' ~ 

010 e 1—11 sao vetores como os procu- 



9.2 Perturbagoes 

Estudaremos agora uma maneira de limitar o efeito de uma pertubagáo nos dados do sistema 
linear, Ax = b, sobre a sua solugáo. 



Teorema 9.1 (da pertubagáo) Se os dados do sistema Ax = b forem pertuhados, isto é, alter- 
ados de uma matriz 5A de um vetor 6b, a resposta será pertubada por um 6x, isto é {A + 6A){x + 
6x) = {b + 6b), tal que, se \\6A\\ \\A-^\\ < 1, entáo 



\\6x 



< 



cond{A) 



(\m\ , m\\\ 

x\\ - l-cond{A) \\6A\\/\\A\\\\\b\\ \\A\\ ) ' 



Demonstragao. 

{A + 6A){x + 6x) = Ax + A6x + 6Ax + 6A6x 
6x = A-^{6b — 6Ax — 6A6x). Assim, 



115x11 < \\A 



-ii 



h + 6b, portanto A6x 

||M||||x|| + ||M||||fe||) . 



— 6Ax — 6A6x, e 
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Como também ||6|| < ||A||||x||, temos que 

\\6x 



{l-\\A-'\\\\6A\\ 



(\m 



\x\ 



< |U-^|| (^^ + ||M| 



Usando a hipótese ||(5A||||A~^|| < 1, temos o teorema, QED. 

Exemplo 2: 

Considere o sistema {B + 6A)x = {b + 6b), onde a matriz de coeficientes, e o vetor de termos 
independentes correspondem á matriz binomial de dimensáo 3, ^B, e ^b \ ^Bl = ^b. Os elementos 
da matriz e do vetor de pertubagáo, 6Al e 6bi, sáo aleatórios. A distribuigáo de cada um destes 
elementos de perturbagáo é independente e é uma fungáo de dois parametros {a,p): Tomemos 
cada elemento da perturbagáo é no cojunto {0, —a, a} com probabilidade, respectivamente, [1 — 
p,p/2,p/2]. 

Com os dados do Exemplo 1, determine um limite máximo para < o; < alfamax que 
garanta a hipótese do teorema da pertubagáo. Faga uma experiéncia comparando hmites de erro 
e pertubagáo da solugáo do sistema proposto. 

Na norma 1, \\6A\\ < 3a, e do Exemplo 1 sabemos que ||-B~^|| < 7. Assim, a condigáo 
||(5A||||A~-'^|| < 1 está garantida para a < 0.04 < 1/21. 

Tomando como experimento de perturbagáo. 



6A 



0.01 -0.01 
0.01 

0.01 -0.01 



e 6b 





-0.01 





a solugao do sistema pertubado é x = [1.0408,0.9388,1.0622]', i.e. 6x 
\\6x\\ =0.16. 

Por outro lado, o Teorema da pertubagáo nos dá o hmite 



[0.04,-0.06,0.06]' e 



\\6x\ 



< 



49 



,0.01 0.03, 



0.35 



\x\ 



1-49*0.03/7' 7 7 ' 

um hmite superior de acordo com o resuhado obtido no experimento. 



9.3 Erro na Fatoracáo LU 



Consideremos a resolugao de um sistema, Ax = b, pelo método de Doohttle, isto é, a decomposigao 
A = LU, e a solugáo do sistema Ly = b e Ux = b. 
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Teorema 9.2 (Wilkinson) A solugáo, afetada pelos erros de arredondamento , pode ser escrita 
como a solugáo exata, x, mais um termo de pertubagáo 6x. Assumiremos que um produto escalar 
em dupla precisáo é afetado de um erro da ordem do erro de passagem para precisáo simples, 
conforme a observagáo 8.1. Nestas condigóes,, a solugáo calculada, {x + 5x), é solugáo exata de 
um sistema. {A + 6A){x + 6x) = b, tal que \\6A\\oo < (2n + l)gu, onde g = maxjj \U¡\. 

Demonstragáo. 

Consideremos as linhas da matriz A já ordenadas de modo a náo haver necessidade de piv- 
oteamentos. A decomposigáo da matriz A será dada por, para i = 1 . . .n 

UU = Ai 
LU' = A' 

ou, para i = 1 . . .n, para 



j=t...n UÍ = A¡-J2MtUi 

k=l 

j=t + l...n M] = {A]-Y.M';Ul)/Ul 

k=l 

Na reahdade, obteremos as matrizes afetadas de erro, para i = 1 . . .n, para 

j=z...n ÜÍ = fl{fld{A¡-Y.M^Üi)) 

k=l 

j = i + l...n M] = fl{fld{{A] - Y: M^ÜD/ÜÍ))) 



fc=i 

Como supomos desprezíveis os termos de 0{ud) (vide observagáo 8.1) temos, para i = 1 . . .n, 
para 

j=t...n ÜÍ = a + SÍ){A>-YM^Ui) 

k=l 

j = t + l...n MÍ = Í1 + S]){A]-Y: M^Ül)/Ü¡ 

k=l 

Portanto, para i = 1 . . .n, para 

j=t...n Y MfÜi + 1UI = UÜ' =A\ + ÜÍ61/{1 + 6t) 

k=l 
i~l 

j=i + l...n Y. M-ÜÍ + MÍÜÍ = Ífi' = A] + M]ÜÍ6]/{1 + 6]) 



k=l 
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Notemos agora que: |(5/|/(1 + S¡) ~ \S¡\ < u, que |M/| < 1, pelo pivoteamento parcial, e 
definindo g = maxjj |f//|, temos que LU = A + E, onde \E¡\ < gu, donde ||-E||oo < ngu. 

Na solugáo do sistema Ly = b e Ux = y calculamos, para i = 1 . . .n. 



e novamente, para i = 1 . . .n. 



y. = fl{fld{{h-Y.L¡y,)/L¡)) 



X. = fl{fld{{y. - Y: U¡xj)/Ui)) . 
j=i+i 



Supondo desprezíveis os termos de 0{ud) 



ou 



y, = {l + S,){h-Y.Lly,)/^ 
i=i 

n 

X, = {l + S[){y,- Y. UÍ¿,)/ÜÍ 
j=i+i 



J2 UlVj = Liy = bi + L¡yiSi/{l + Si) 
i=i 

n 

^f//x, = Ü,x = y, + UÍxX/{l + S'i) 
j=i 



isto é 

{L + SL)y = b 
{Ü + SU)x = y 

onde para as matrizes diagonais SL e SU, temos \SLl\ < u e \SU¡\ < gu. 
Em suma, 

b = {L + SL)y = {L + SL) {Ü + SU)x 

= {LÜ + LSU + SLÜ + SLSU)x 

= {A + E + LSU + SLÜ + SLSU)x 

= {A + SA)Í 

Desprezando o termo SLSU , de 0{ud), SA = E + LSU + SLU , donde 

I |M| loo = 11^11 + \\LSU + SLUÜ\ I < ngu + {n+ l)gu 
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pois 



L6U + 6LU oo < 




















" 1 


" 




gu 


' 


u 


" 




' 9 


< 










+ 












_ 1 


1 _ 




gu _ 







u _ 




_ 




" gu 





' 


' gu gu' 






< 






+ 










_ gu 


gu 


. 


_ gu _ 








' 2gu 


gu 1 






< 






= (n + l)gu QED. 






_ gu 




2g 


u J 















Para frisar a importáncia do uso conjugado do método de Doolittle e dupla precisáo daremos, 
sem demonstragáo, a versáo do teorema ora demonstrado para o método de Gauss, que equivale 
ao método de Doolittle sem o recurso da dupla precisáo: Definindo h = maxij^k \ ^^¿|, teríamos 
||-E||oo < Oirí^hu) e ||(5yl||oo < 0{n^hu). Estes resultados tornam evidente a obrigatoriedade de 
calcularmos os produtos escalares envolvidos no processo, sempre em dupla precisáo, a menos que 
tratemos de sistemas de pequena dimensáo. 

Além do pivoteamento parcial, isto é, da escolha como elemento pivó do maior elemento em 
módulo na coluna, poderíamos usar o pivoteamento total, isto é, encolher como elemento pivo 
o maior elemento em módulo em qualquer hnha ou coluna utihzável (isto é, poderiamos tomar por 
pivó da transformagáo ^~^A — )> '^A, um elemento em arg maxkKijKn \ ^'^Aj. Usando pivoteamento 
total, podemos demonstrar a existéncia de um hmite superior para a constante h, supondo que 
a matriz original A é tal que |A¿| < 1, da ordem de 0(n*^^/^^'°*^")), contra hmites de 0(2") para 
pivoteamento parcial. O uso do pivoteamento total é todavia, extremamente custoso, pois exige 
a cada passo da O^ri^) comparagóes, contra 0{n) no pivoteamento parciaL Ademais, estes hmites 
de h tendem a ser extremamente pessimistas, principalmente se A for bem equihbrada. 



Exercícios 



1. Estabihdade de Fatoragoes Simétricas. 



(a) Adapte os resuhados de estabihdade da faoragao LU para o caso particular da fatoragao 
de Cholesky. 

(b) Qual o número de condigáo de uma matriz ortogonal? Usando a relagao entre o fa- 
tor triangular da fatoragáo QR e a fatoragáo de Cholesky, discuta a estabihdade da 
fatoragáo QR. 
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2. Prove que || ||i, || II2 e || ||oo sáo efetivamente normas. 

3. Desenhe em R^ a regiao ||a;|| < 1, para as normas 1, 2 e 00. 

4. Prove o Lema 1. 

5. Equivaléncia das normas 1, 2 e 00: Prove em R" que 

T < Ti ^TÍT 

|*^||oo _ 11*^111 — "'ll'^lloo 

T ^ T o ^ Tí / T 

|*^||oo _ ||*^||z _ '*' 11*^1100 

||A||2 = A, onde A^ é o maior autovalor de A'A. 

I 4' /i|lo — II 4I|2 
jí jí 2 — -^ 2 



(aj 
(b) 
(c) 
(d) 
íe) 



6. Calcule computacionalmente cond{A), para as matrizes de teste "T, A = "^B e A = '^H, 
para n = 2, 4, 8, 16. Estime graficamente o crescimento do número de condigáo das matrizes 
teste em fungáo da dimensáo. 

7. Resolva computacionalmente os sistemas teste ^Tx = ^t, ^Bx = ^b e ^Hx = ^h. 

(a) Pelo método de Gauss com pivoteamento parcial; 

(b) Pelo método de Gauss com pivoteamento total; 

(c) Pelo método de Doolittle com pivoteamento parcial e dupla precisáo; 

(d) Pelo método de Doolittle com pivoteamento total e dupla precisáo. 

No item c, verifique se o erro final está dentro do esperado, em fungáo da unidade de erro 
do SPF utilizado. 



Capítulo 10 
MUDANQA de BASE 



Atualizagoes de Posto 1 



Em Otimizagáo, principalmente em Programagáo Linear, o termo base significa uma matriz 
quadrada de posto pleno. Estudaremos agora o problema de mudanga de base, aqui posto na 
seguinte forma: Seja A a nova base, obtida da base original, A, substituindo-se a coluna A^ por 
uma nova coluna, a, isto é A = [A^, . . . A^^^, a, A^^'^, . . . A"]. Se já dispusermos da inversa de A (ou 
na prática mais comumente de uma fatoragáo de A), como atualizar a inversa (ou a fatoragáo)? 
Isto é, como, a partir da inversa de A obter a inversa de A com um trabalho muito menor que 
o necessário para reinverter A, i.e., computar a inversa de novo, sem aproveitar a informagáo 
contida em A~^ ou, na fatoragáo A = LU on A = QR ? 



10.1 Fórmulas de Modificagáo 

Teorema 10.1 (fórmula geral de modificagáo) Dada A, n x n e inverswel, V = A~^ , 6A, 
n X n, e a E R suficientemente pequeno, podemos fazer a expansao 

oo 

{A + a5A)-^ = V + Y.i-a)^ {V 5A)^ V . 

k=l 

Demonstragáo: 

Em virtude da regra de Cramer podemos, para a G [0, alphamax], escrever a série de Taylor 
para cada elemento de {A + a5A)~^, 

oo k 

{A + 5A)-' = V + Y.iy^^ 
k=i ^- 

Para determinar a matriz que dá a perturbagáo de ordem k da inversa 5^V observemos que 

2 3 

(A + a5A){V + a5V + — 5^V + ^5'^V + ...) = ! 
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donde 

a{6A V + A6V) + a^^A 6^V + 6A 6V) + «^(^A íV + ]-6A 6^V) + . . . = 

2 6 2 

de modo que 

6V = -V 6AV , (5V = 2V 6AV 6AV , 6^V = {V 6A)^ V ,... 
6^V = (-1)^ k\ iy 6A)^ V . 

Substituindo esta fórmula geral de 6^V na série de Taylor de {A + a 6A), segue diretamente o 
teorema. QED. 

Teorema 10.2 (fórmula de Sherman e Morrison) Dada A, n x n e inverswel, u e w n x 1, 
e a & R* , entáo a inversa de A = A + auw' é dada por 

Á~^ = A-^ + (3A-\w'A-^ , 

/3 = -{a-^+w'A-\)-^ . 

Demonstracáo: 

Da fórmula geral de modificagáo, supondo que a < alfamax, e das propriedades do operador 
trago (exercício 1.5) temos: 



(A 


+ auw' 


r' = 








= v + 


co 

> ]{-a)^{yuw')V 

k=l 






= v - 


aVuw'V)_^{- 

h 


-atr{Vuw 


''))' 




= v - 


aVuw'V 


) 






1 + atiiVuw' 






= v - 


Vuw'V 







a ^ + ti^Vuw') 
= V - {a-^ + w'Vu)-^ Vuw'V . 

o que prova o teorema para < o; < alphamax. 

A fórmula de Sherman e Morrison é todavia uma identidade que podemos provar diretamente 
para qualquer o; > | {a^^ + w'A^^u) ^ 0: Usando a formula de Sherman e Morisson para 
desenvolver a identidade Ar^A = I, obtemos. 



{V + (5Vuw'V){A + auw') = I ^ 

/3Vuw' + a/3{Vuw')'^ = -aVuw' <^ 
(q;~^ + w'Vu)aVuw' = -f^-^aVuw' 
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sendo esta última identidade trivialmente verdadeira. Q.E.D. 
Exemplo 1: 
Tomando 



A 



1 

2 1 



u 



w 



3 

4 



, a = 1 , 



podemos calcular a inversa de 



A 



"10" 


+ 


"34" 




"44" 


2 1 




_ 6 8 _ 




_ 8 9 _ 



calculando 



Á-' 



1 
-2 1 



/3 = -(1 + 3)-^ = -1/4 
- (1/4) 



-5 4 




9/4 -1 
-2 1 



Observagáo 10.1 Note que uma mudanga de base poderia ser feita pela fórmula de Sherman e 
Morrison, pois se A = [A^, . . . A^~^, a, A^~^^, . . . A"], entáo, A = A + {a — A^) Ij, de modo que a 
inversa da nova base seria, tomando V = A~^, 

i-i = V + /3V{a-A^)IjV 
= V + V/3{a-A^)Vj 

onde /3 = -{l- IjV{a - A^))-'^ = -{Vja)-^ . 



Exemplo 2: 

Se 



A 



1 

2 1 



1 
3 



e J = 2 , 



entáo a inversa da nova base, A ^, pose ser calculada como segue 



/5 



-2 1 



1 
3 



v-l 



i-1 



1 1 

2 3 

1 

-2 1 



1 

-2 



" 1 ■ 


+ 


" -1 " 


-2 1 








-1)( 



-2 1 



1 
3 





1 

-2 1 



-2 1 
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10.2 Atualizacoes Estáveis 

O Algoritmo de Bartels e Golub, que apresentaremos a seguir, nos dá uma atualizagáo estável 
da fatoracáo LU de uma matriz. Tomando T = L~^, A = LU, 

TÁ=[U\... U'-\ Ta, W+\ ...U'']. 

Consideremos agora a permutagáo de colunas que comuta a s-ésima e a n-ésima colunas, Q 
Esta permutagáo leva A em A = AQ, e TA em 



TAQ 



U¡ ... Ur Tia Ui 



s+l 



H 







, 

ul ... 

••. 






... 



Ulzl T,^,a U!+¡ 
Tsa U¡+'^ 
T,+^a U¡Xl 



T„a 



U{' 

Tjn 
^ s-l 

jjn 

^ s 
Tjn 
^ s+l 

jjn 



Q 



f/r^ ui 



s+l 



rn-l 



jjn-, jjn j.^^ 



U 



s-1 
s-1 












JJS + 1 
JJS + 1 

^ s 

JTS + 1 
^S + l 








rrn-l 
^s-l 



U. 



s-1 

í/r 



Ts-ia 
Tm 



U^+i f^s'+i Ts+ia 



Un-l ^ji-l Tn-ia 







jjn 
^ n 



T„a 



H é uma matriz de Hessenberg superior, isto é, apenas os elementos paralelos á diagonal 
principal podem ser náo nulos no triángulo inferior. O método de Bartels e Golub consiste na 
aplicagáo de método de Gauss, com pivotamento parcial, á matriz H. 

Observemos que a aplicagáo do método de Gauss á matriz H é extremamente simples, pois 

1. As primeiras s — 1 etapas de transformagáo náo sáo necessárias, pois nas colunas 1, ... s — 1, 
if já é triangular superior. 

2. As etapas j = s, . . .n — 1, que transformam 



^H = '-^H -^ 'H ^ ... ^H ^ ... ^-^H = Ú 
resumem-se em 

(a) Permutar as linhas j e j + 1 se | ^^^Hj^i \ > \ ■'^^Hj \, obtendo uma nova matriz ^^^H. 
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(b) Calcular um único multiplicador "'A^j+i = "' ^-^j+i/ "' ^-^j- 

(c) Atualizar uma única linha ■'Hj^i = ^^^Hjj^i — ■' Njj^-^ ^^^Hj. 



A aplicagao do metodo de Gauss nos dá a fatoragao H = LU, onde H = RH é a matriz obtida 
de H pelas permutagóes de linhas reahzadas durante a triangularizagáo. Assim, 



H 

Á 

V = Á~^ 



RH = RTÁ = LU , donde 

LR'LÚ , ou 

Ú^^fRT 



Após uma seqüéncia de mudangas de base, nossa representagáo da inversa teria a forma 



V= ^U'^^T^R...^T^RT 



10.3 Preservando Esparsidade 

O algoritmo de Saunders, que agora examinamos, pode ser visto como uma adaptagáo do 
algoritmo de Bartels e Golub visando aproveitar a prévia estruturagáo da base pelo algoritmo P4. 

Suponhamos termos a fatoragáo A = LU obtida pela aphcagáo do método de Gauss (sem 
pivotamento de hnhas mas com possíveis permutagoes de colunas) á matriz A previamente estru- 
turada pelo P4, por exemplo 

Exemplo 3: 





1 


2 


3 


4 


5 


6 


7 8 


9 


10 


1 


X 








2 




X 




X 












3 


X 


X 


X 















4 






X 


X 












5 


X 


X 


X 




6 


X 




X 







X 






X 


7 




X 




X 


X 




X X 




X 


8 




X 











X X 




X 


9 




X 




X 




X 


X 


X 


X 


10 


X 








X 




X 


X 






No Exemplo 3, "x" indica as posigoes originalmente náo nulas de A e "0" indica as posigoes 
preenchidas durante a triangularizagáo. 
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Seja f/ a matriz obtida de U pela permutagao simétrica, Q'UQ, que leva os espinhos, preser- 
vando sua ordem de posicionamento, para as últimas colunas á direita. No Exemplo 3, teríamos. 









1 


2 


3 


5 


6 7 9 


4 8 


10 




1 


X 


















2 




X 








X 








3 






X 






X 








5 








X 








ü = 


= Q'UQ = 


6 
7 
9 










X 

X 

X 


X 


X 

X 




4 












X 








8 












X 


X 






10 














X 



Esta matriz tem estrutura U 



D E 
F 



, onde D é diagonal e F é triangular superior. 



O algoritmo de Saunders atualiza a decomposigáo A = LQUQ' para a nova base 



A 
ÁQ 



A^ ... A'-^ a A'+^ ... A" 
A^ ... A'-^ A'+^ ... A" a 



ou 



como segue: 



1. Forma, pela permutagao UQ, e substituigao da última coluna, 

W = TÁQ = Q'TAQQ =\Ü^ ... Ü'-^ Ü'+^ ... ¿/" Q'Ta 



2. Forma, pela permutagáo simétrica. 



W = Q'TÁQ = 
{Wi)' ... {Ws-i)' {Ws+i)' 



= Q'Q'TÁQQ = W 

■ ■ ■ {Wn)' {Ws)' " 



Note que W tem uma estrutura do tipo 



W 



D E W-^_,_, 

" " n—c:n—l 

w, w:: 



onde D é diagonal e W^ só pode ter elementos nao nulos nas últimas posigoes á direita (sob 
Fe W^" ). 
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3. Triangulariza W, isto é, N, pelo método de Gauss com pivotamento parcial. A matriz Á^, 
constituída de F e dos últimos elementos da linha Wg e da coluna W"', é denominada núcleo 
de W. Nas rotinas numéricas pode ser conveniente guardar apenas o núcleo, c x c, c « n 
como uma matriz densa na memória principal, enquanto o resto da matriz, usada apenas 
para leitura, pode ser guardada em uma representagáo esparsa, e na memória secundária. 



No Exemplo 3 teríamos, ao mudar a coluna 3 da base A, 



A 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


1 


X 




y 
















2 




X 




X 














3 


X 


X 


















4 






y 


X 














5 


X 






X 


X 












6 


X 










X 








X 


7 




X 


y 


X 


X 




X 


X 




X 


8 




X 










X 


X 




X 


9 




X 




X 




X 




X 


X 


X 


10 


X 




y 




X 






X 


X 





w 





1 


2 


5 


6 


7 


9 


4 


8 


10 


3 


1 


X 


















y 


2 




X 










X 








3 














X 








5 






X 
















6 








X 










X 




7 










X 






X 


X 


y 


9 












X 






X 




4 














X 






y 


8 
















X 


X 




10 


















X 


y 
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W 





1 


2 


5 


6 


7 9 


4 


8 


10 


3 


1 


X 
















y 


2 




X 








X 








5 






X 














6 








X 








X 




7 










X 




X 


X 


y 


9 










X 






X 




4 












X 






y 


8 














X 


X 




10 
















X 


y 


3 












X 










Após a triangularizagao de W teremos, sendo R as permutagoes de linhas realizadas ao trian- 
mlarizar W, 



RW = LU = RQ'Q'TAQQ portanto 

Á = LQQR'LÚQ'Q' e 
i-i = QQÚ-^TRQ'Q'T . 

Em geral, após uma seqüéncia de mudangas de base, nossa representagáo da inversa terá a 
forma 

M"^ = Q^Q... ^Q ^f/-i '^T ^R ^Q' ... ^T^R^QQT . 



Observagáo 10.2 



1. A regiao dos preenchimentos em ^W está restrita ao triángulo superior do núcleo e sua 
última hnha. 

2. A matriz ^L só terá elementos nulos na úhima hnha do núcleo. 

3. A cada nova mudanga de base, a dimensáo do núcleo: 

(a) Aumenta de 1, se a coluna que sai da base é uma coluna triangular da base original. 

(b) Permanece constante, se a coluna que sai da base é um espinho, ou uma coluna já 
anteriormente substituída. 

4. Cada mudanga de base aumenta um termo na fatoragáo da base, o que leva a uma gradativa 
perda de eñciéncia e acúmulo de erros. Depois de um número pré-determinado de mudangas 
sobre uma base original, ou quando o acúmulo de erros assim o exigir, devemos fazer uma 
reiversáo, i.é., reiniciar o processo aphcando a P4 e triangularizando a próxima base desejada. 
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10.4 Preservando Estrutura 

Consideraremos agora o problema de atualizar a fatoracáo A = QU de uma base na forma angular 
blocada. Conforme argumentamos no final do capítulo 7, estamos apenas interessados no fator 
U, sendo as transformagóes ortogonais descartadas logo após o seu uso. No que segue, a coluna a 
sair da base será a coluna outj do bloco outk, 1 < outk < h + 1. Em seu lugar será introduzida 
na base uma coluna com a estrutura do bloco ink, 1 < ink < h + 1. 

Apresentamos agora um procedimento de atualizagáo de U por blocos, bup{), que usa explici- 
tamente a estrutura angular blocada da base A e do fator U. Este procedimento será descrito em 
termos das operagoes simples em cada bloco apresentadas na segáo 7.5. 

Em bupO consideraremos cinco casos distintos: 

Caso I ink ^ outk, ink ^ h+ 1, outk ^ h+ 1. 
Caso II ink = outk, ink ^ h + 1. 
Caso III ink ^ h + 1, outk = h + 1. 
Caso IV ink = h + 1, outk ^ h + 1. 
Caso V ink = outk, ink = h + 1. 

Vejamos em detalhe o caso I, quando ink and outk sáo blocos diagonais distintos, como 
mostrado na figura 1. Neste caso os únicos ENN's na coluna saindo estáo no bloco o™^iJ, J. 



Analogamente os únicos ENN's na coluna entrando na base, a, estao em 



ink 



a. 



Definimos y = A'a, e u = Q'a = U ^A'a = U *í/. Notamos que os vetores y e u preservam a 
estrutura blocada da base. Assim, os ENN's em u estáo nos blocos *"^m e ^~^^u. 



ink 
h+1 



U 



U 



ink'í/' ink'íx/' 

s 



ink„ 
h+l 



y 



Para atuaUzar U removemos a coluna outj do bloco outk, e inserimos u como a última col- 
una de ^'^f/. Depois apenas temos que reduzir U a uma matriz triangular superior através de 
transformagóes ortogonais. 



1. Leve 

2. Leve 



outk '\r outk Tir 
^+^M S 



de Hessenberg a triangular superior. 

a triangular. 

3. Insira a primeira linha de ^~^^U , como a última hnha de *"'=f/. 

4. Insira a última linha de ""^'^f/ como a primeira de ^+^f/. 
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5. Leve S de Hessenberg a triangular superior. 

Os outros casos sáo bastante similares. Os esquemas nas figuras 2, 3 e 4 sáo os análogos do 
esquema na figura 1, e dáo uma descrigáo sumária de cada um dos casos. 

Estes sáo os passos para o caso I : 



1. No nó ink, compute y 



ink 



í ninkyt „ink „ n.fe+l 



{C 



ink\t„ink 



pTime = m{ink)n{ink) + m{ink)n{b + 1) < ^dhmax^ , INC 
2. No nó ink, compute a transformacáo inversa parcial 



. 



u 



ink 



'í/'ink \x/'ink 

/ 



-I -t r 



,ink 
6+1 



Entao insira ■u*"'^ como a última coluna de 1^*"-^. 

pTime = {l/2)n{ink)^ + n{ink)n{b + 1) < {3/2)dbmax^ , INC = . 

3. Do nó ink envie z ao nó 0. 

pTime = , INC = n{b + 1) < dbmax . 



4. No nó compute u^^^ = S *z. 

pTime = (l/2)n(6+ 1)2 < {l/^^dbmax"^ , INC 



0. 



5. (a) No nó outk, remova a coluna V,""*^- de \/°"*'=. Entao reduza 



senberg para triangular superior. 



V 



outk 



w 



outk 



de Hes- 



(b) No nó 0, reduza 



u 



6+1 



S 



para triangular superior. 

Observe que as operagóes nos passos 5a e 56 sáo independentes, portanto 
pTime = ^n^ink)"^ + An{ink)n{b + 1) A 2n{b + 1)^ < Gdbmax'^ , INC = 0. 

6. Do nó envie o vetor 5*1^, ao nó ink, onde ele é inserido como a última coluna de W^*"^. Do 
nó envie o elemento u\'^^ ao nó ink, onde ele é inserido como U^n(ink)+iMink)+i- -Do nó outk 
envie vetor W^/^'^^^-. , ao nó 0, onde ele é inserido como a primeira linha de S. 

pTime = , INC = 2n{b + 1) + n{outk) < Sdbmax. 

7. No nó 0, reduza S de Hessenberg para triangular superior. 
pTime = 2n{b + 1)^ < 2dbmax^ , INC = 0. 

Estes sáo os passos para o caso II : 

Os Passos 1 — 5 sáo exatamente como no caso I. 



6. (a) Do nó ink envie V;Y^fc),n(mfc) e W^^"f„fc),. ao nó 0. 
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(b) No nó reduza para triangular superior a matriz 2 x n{b 



'írink jx/'ink 

n(ink) ,n{ink) n{ink),* 



"1 



-51,. 



pTime = 4n{b + 1) < Adbmax , INC = n{b + 1) < dbmax. 

7. (a) Do nó envie o vetor modificado [ KÍSfc),„(mfc) ^Ífnk),, 
(b) No nó 0, reduza S de Hessenberg para triangular superior. 
pTime = 2n{b + 1)^ < ^dbmax"^ , INC = n{b + 1) < dbmax. 



de volta ao nó ink. 



Estes sao os passos do caso III : 

Os passos 1 — 4 sáo exatamente como no caso I. 



5. (a) No nó k = 1 : b remova a coluna W^^^^^j de W'^. 



(b) No nó reduza 



u 



b+l 



S 



para triangular superior. Remova S,^outj de S. 



pTime = 2n{b + 1)^ < 2dbmax^ , INC = 0. 

6. Do nó envie ao nó ink, u\'^^ para ser inserido em V'"'^'' como Kí{¿nfc)+i,n(mA:)+i) ^ '^'i,, para 
ser inserido como a última coluna de W^"'''. 

pTime = , INC = n{b + 1) < dbmax. 

7. No nó 0, reduza S de Hessenberg para triangular superior. 
pTime = 2n{b + 1)^ < 2dbmax^ , INC = 0. 



Estes sao os passos do caso IV : 



1. No nó k=l:b, compute y'' = {B'^Ya'' e x'' = {C'^Ya'^. 

pTime = /\\ m{k)n{ink) + m{k)n{b + 1) < ^dbmax"^ , INC = 0. 

2. No nó k=l:b, compute a transformagáo inversa parcial 



u 



yk ^rk 

/ 



,ink 



X 



e insira u'' como a última coluna de W''. 

pTime = A\ (l/2)n(A;)2 + n{k)n{b + 1) < {S/^^dbmax"^ , INC = 0. 

3. Do nó k=l:b envie z'' ao nó 0, onde acumulamos z = J2\ z'^. 
pTime = b n{b + 1) < fe dbmax , INC = b n{b + 1) < 6 dbmax. 

4. No nó compute u^^^ = S^'^z, e insira -u*^^ como a última coluna de S. 
pTime = {l/2)n{b + 1)^ < {l/^^dbmax"^ , INC = 0. 
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5. remova a coluna V,°o*^- de V"'^*'' , e reduza 



outk Txroutk 



youtk ^r 



para triangular superior. 



pTime = 2n{outkY + 4:n{outk)n{b + 1) < Qdhmax^ , INC = 0. 



6. Envie o vetor W°y^^¡^^ , do nó outk ao nó 0, onde o inserimos como a primeira linha de S, e 
reduza S a triangular. 
pTime = 2n{b + 1)^ < 2dhmax^ , INC = n{h + 1) < dhmax. 



Estes sao os passos do caso V : 

Os passos 1 — 4 sáo exatamente como no caso IV. 



5. No nó k=l:b, remova a coluna W^^^^, de W°^^^ , e insira u^ como a última coluna de W^ 
No nó remova a coluna S,^outj de S* , e insira -u^+^ como a última coluna de S. 
pTime = , INC = 0. 



6. No nó 0, reduza S de Hessenberg para triangular superior. 
pTime = 2n{h + 1)^ < 2dhmax^ , INC = 0. 



A complexidade do procedimento de atualizagao de U por blocos, hup(), é dada pelo seguinte 
teorema: Como na segáo 7.4, dhmax = max.'¿^-¡^n{k). 



Teorema 10.3 A complexidade de hup{), desconsideramos termos de ordem inferior, tem os 
seguintes limitantes superiores para tempo de processamento e comunicagáo: 



Caso 


Processamento 


Comunicagáo 


I 


I2dhmax'^ 


Sdhmax 


II 


12dhmax'^ 


3dhmax 


III 


Sdhmax"^ 


2dhmax 


IV 


12dhmax'^ + h dhmax 


h dhmax 


V 


Gdhmax"^ + h dhmax 


h dhmax 



Caso o ambiente permita comunicagoes em paralelo, podemos substituir h por log{h) nas 
expressóes de complexidade. 
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1 1 1 
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1 1 








C' 


outk = 1 
1 1 1 
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1 1 

ink = 2 
1 1 




c^ - 






1 
B^ 


c^ 
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1 1 







-^ 


V' 








u 




— 


w^ 




outj 


= 2 






V 


'2 




w 








v^ 




w^ 








s 













u 
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Caso IV 



Caso V 
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Histórico 

Matlab, de Matrix Laboratory, é um ambiente interativo para computacáo envolvendo matrizes. 
Matlab foi desenvolvido no inicio da década de 80 por Cleve Moler, no Departamento de Ciéncia 
da Computagáo da Universidade do Novo México, EUA. 

Matlab coloca á disposigáo do usuário, num ambiente interativo, as bibliotecas desenvolvidas 
nos projetos LINPACK e EISPACK. Estes projetos elaboraram bibliotecas de domínio público 
para Álgebra Linear. LINPACK tem rotinas para solugáo de sistemas de equagoes lineares, e 
EISPACK tem rotinas para cálculo de autovalores. Os manuais destes projetos sáo portanto 
documentagáo complementar á documentagáo do Matlab. 

Versóes posteriores de Matlab, atualmente na versáo 4.0, foram desenvolvidas na firma comer- 
cial MathWorks Inc, que detém os direitos autorais destas implementagóes. As versóes recentes do 
produto Matlab melhoram significativamente o ambiente interativo, incluindo facilidades gráficas 
de visualizagáo e impressáo; todavia a "Linguagem Matlab" manteve-se quase inalterada. Existem 
vários interpretadores da linguagem Matlab em domínio publico, como Matlab 1.0, Octave e rlab. 
Existem também outros interpretadores comerciais de Matlab, como CLAM. Existem ainda várias 
Tool Boxes, bibliotecas vendidas pela MathWorks e por terceiros, com rotinas em Matlab para 
áreas específicas. 

Usaremos a grafia de nome próprio, Matlab, como referéncia a linguagem, o nome em maiúsculas, 
MATLAB, como referéncia ao produto comercial da MathWorks, e a grafia em minúsculas, mat- 
lab, como referéncia a um interpretador genérico da linguagem Matlab. 

O Ambiente 

Para entrar no ambiente matlab, simplesmente digite "matlab". O prompt do matlab é >> , que 
espera por comandos. Para sair use o comando quit. Dentro do ambiente matlab, um comando 
precedido do bang, !, é executado pelo sistema operacional, assim: usando !dir ou !ls ficaremos 
sabendo os arquivos no diretório corrente, e usando !edit , !vi ou !emacs , seremos capazes de 
editar um arquivo. Normalmente Matlab distingue maiúsculas de minúsculas. 
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O comando help exibe todos os comandos e símbolos sintáticos disponíveis. O comando 
help nome com fornece informagoes sobreo comando de nome nomecom. O comando diary nomearq 
abre um arquivo de nome nomearq, e ecoa tudo que aparece na tela para este arquivo. Repetindo 
o comando diary fechamos este arquivo. O formato dos números na tela pode ser alterado com 
o comando f ormat. 

Os comandos who e whos Ustam as variáveis em existéncia no espago de trabalho. O comando 
clear Umpa o espago de trabalho, extinguindo todas as variáveis. O comando save nomearq 
guarda o espago de trabalho no arquivo nomearq, e o comando load nomearq carrega um espago 
de trabalho previamente guardado com o comando save. 

Em Matlab há dois terminadores de comando: a vírgula, , , e o ponto-e-vírgula, ; . O resultado 
de um comando terminado por vírgula é ecoado para a tela. O terminador ponto-e-vírgula náo 
causa eco. Resultados ecoados na tela sáo atribuídos a variável do sistema ans (de answer, ou 
resposta). O terminador vírgula pode ser suprimido no último comando da Unha. Para continuar 
um comando na hnha seguinte devemos terminar a hnha corrente com trés pontos, . . . , o símbolo 
de continuagáo. O sinal de porcento, /d , indica que o resto da hnha é comentário. 

Matrizes 

O tipo básico do Matlab é uma matriz de números complexos. Uma matriz real é uma matriz que 
tem a parte imaginária de todos os seus elementos nula. Um vetor hnha é uma matriz 1 x n, um 
vetor coluna uma matriz n x 1, e um escalar uma matriz 1x1. 

As atribuigóes 
A = [1, 2, 3; 4, 5, 6; 7, 8, 9] ; ou 
A = [12 3; 4 5 6; 7 8 9];, 
sáo equivalentes, e atribuem á variável A o valor 

1 
A 



2 

5 



Matrizes sáo deUmitadas por colchetes, elementos de uma mesma hnha sáo separados por 
vírgulas (ou apenas por espagos em branco), e UnUas sáo separadas por ponto-e-vírgula (ou pelo 
caracter nova-UnUa). O apóstrofe, ' , transpóem uma matriz. É fácil em Matlab compor matrizes 
blocadas, desde que os blocos tenUam dimensóes compatíveis! Por exemplo: 



A 



"12" 
3 4 


, B = 


' 5 ■ 
6 


, c = 


' 7' 
8 










9 



D=[A,B;C']; D 



1 
3 

7 



Cuidado ao concatenar matrizes com os espacos em branco, pois estes sao equivalentes a vírgulas, 
separando elementos. Assim: [1,2+3]==[1 5] mas [1,2 +3] ==[1,2, 3] . 
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Há várias fungóes para gerar matrizes e vetores especiais: zeros (m , n) , ones (m , n) e rand (m , n) 
sáo matrizes de dimensáo m x n com, respectivamente, zeros, uns, e números aleatórios em [0,1]. 
O vetor i : k : j é o vetor linha [i,i + k,i + 2k, . . .i + nk] , onde n = max m \ i + mk < j. Podemos 
suprimir o "passo" k = 1, escrevendo o vetor i:l:j simplesmente como i:j. Se f é um vetor, 
diag(v) é a matriz diagonal com diagonal v, se A é uma matriz quadrada, diag(A) é o vetor com 
os elementos da diagonal principal de A. A matriz identidade de ordem n é eye(n), o mesmo que 
diag(ones(l,n)). 

A(i, j) é o elemento na i-ésima linha, j-ésima coluna de A, m x n. Se vi e vj sáo vetores de 
índices em A, i.e. vetores hnha com elementos inteiros positivos, em vi náo excedendo m, e em vj 
náo excedendo n, A(vi,vj) é a sub-matriz do elementos de A com índice de hnha em vi e índice 
de coluna em vj. Em particular A(l:m,j), ou A(: ,j) é a j-ésima coluna de de A, e A(i,l: j), 
ou A(i, : ), é a i-ésima hnha de A. Exemplo: 



A 



11 12 13 " 




" 3 ■ 
1 




" 2 " 
3 


21 22 23 


vi = 


vj = 


31 32 33 







A{vi,vj] 



32 33 

12 13 



As operagóes de adigáo, subtragáo, produto, divisáo e poténcia, + -*/", sáo as usuais da 
álgebra hnear. O operador de divisáo á esquerda, \, fornece em x = A\b ; uma solugáo x \ A*x = b. 
O operador de divisáo a direita, / , fornece em a = b/A ; uma solugáo x \ x * A = b. Quando o 
sistema é bem determinado isto é o mesmo que, respectivamente, inv(A) *b ou b*inv(A) . (Quando 
o sistema é super-determinado x é uma solugáo no sentido dos mínimos quadrados.) Os operadores 
aritméticos de produto, poténcia e divisao tem a versáo pontual, . * . " . / . \ , que sao executados 
elemento a elemento. Exemplo: 



X 



5 5 5 



y 



-1 1 



X . * y 



-5 5 



Matlab é uma hnguagem declarativa (em oposigáo a imperativa) onde as variáveis sáo declaradas, 
dimensionadas e redimensionadas dinamicamente pelo interpretador. Assim, se A presentemente 
náo existe, A=ll; declara e iniciahza uma matriz real 1x1. Em seguida, o comando A(2,3)=23; 
redimensionaria A como a matriz 2 X 3 [11, 0, 0; 0, 0, 23]. A nome da matriz nula é [ ]. A 
atribuigáo A( : ,2) = [] ; anularia a segunda coluna de A, tornado-a a matriz 2x 2 [11, 0; 0, 23]. 



Controle de Fluxo 



Os operadores relacionais da hnguagem sao <<=>>= == ~=, que retornam valor ou 1 conforme 
a condigáo seja verdadeira ou falsa. Os operadores lógicos, nao e cm, sáo, respectivamente, ~ & I . 

Matlab possui os comandos de fluxo for - end, while - end e if - elseif - else - end, que tem a 
mesma sintaxe do Fortran, exemphflcada a seguir. Lembre-se de náo escrever a palavra elseif 
como duas palavras separadas. 
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for i=l:n if(x<0) fatn=l; 

for j=l:m signx=-l; while(n>l) 

H(i,j)=l/(i+j-l) ; elseif(x>0) fatn=fatn*n; 

end signx=l; n=n-l; 

end else end 

signx=0; 
end 

Uma consideragáo sobre eficiéncia: Matlab é uma linguagem interpretada, e tempo de inter- 
pretagáo de um comando simples pode ser bem maior que seu tempo de execucáo. Para tornar 
o programa rápido, tente operar sobre matrizes ou blocos, evitando loops explícitos que operem 
numa matriz elemento por elemento. Em outras palavras, tente evitar loops que repetem um 
comando que atribui valores a elementos, por atribuigóes a vetores ou matrizes. As facilidades no 
uso de vetores de índices, os operadores aritméticos pontuais, e fungoes como max , min , sort , 
etc, tornam esta tarefa fácil, e os programas bem mais curtos e legíveis. 

Scripts e Funcoes 

O fluxo do programa também é desviado pela invocagáo de subrotinas. A subrotina de nome 
nomsubrdeve estar guardada no arquivo nomsubr .m ; por isto subrotinas sáo também denominadas 
M-arquivos (M-files). Há dois tipos de subrotinas: Scripts e Fungóes. 

Um script é simplesmente uma seqüéncia de comandos, que seráo executados como se fossem 
digitados ao prompt do matlab. Subrotinas podem invocar outras subrotinas, inclusive recursiva- 
mente. 

Um M-arquivo de fungáo em Matlab comega com a declaragáo da forma 
[psl, ps2, . . . psm] = nomefunc( pel, pe2, . . . pen ) 

A lista entre parénteses é a lista de parámetros de entrada da fungáo, e a lista entre colchetes sáo 
os parámetros de saída. Parámetros sáo variáveis locais, assim como sáo variáveis locais todas as 
variáveis no corpo da fungáo. 

Ao invocarmos a fungáo nomefunc com o comando 
[asl, . . . asm] = nomefunc(ael, . . . aen) ; 

Os argumentos de entrada, ael, . . . aen, sáo passados por valor aos (i.e. copiados nos) parámetros 
de entrada, e ao fim do M-arquivo, ou ao encontrar o comando return , os parámetros de saída 
sáo passados aos argumentos de saída. Exemplos: 

function [mabel, mabin] = vmax(v) 

% procura o maior elemento, em valor absoluto, 

y„ dentro de um vetor, linha ou coluna. 

yoObs: Esta funcao NAO segue os conselhos 
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"/opara operar sobre blocos, e nao elementos! 

[m,n]=size(v) ; 

if( m ~= 1 & n ~= 1 ) 

erro; 
else 

K=max ( [m , n] ) ; 

mabel= abs(v(l)); mabin= 1; 

for k=2:K 

if ( abs(v(k)) > mabel ) 

mabel= abs(v(k)); mabin= k; 
end'/oif 
end7of or 
endZelse 

Para referir-mo-nos, dentro de uma fungáo, a uma variável externa, esta deve ter sido declarada 
uma variável global com o comando global. A forma de especificar uma variável como global muda 
um pouco de interpretador para interpretador, e mesmo de versáo para versáo: Diga help global 
para saber os detalhes de como funciona a sua implementacáo. 
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